Импорт phyloseq объекта и библиотек
Датасет - хроносерия разложения соломы - от фактор Day 10ти уровней
library(phyloseq)
library(tidyverse)
library(ggpubr)
library(ampvis2)
library(ANCOMBC)
library(heatmaply)
library(compositions)
library(igraph)
library(WGCNA)
require(DESeq2)
require(phyloseq)
ps.f <- readRDS("psf2")
#phyloseq object to ampvis2 object
#https://gist.github.com/KasperSkytte/8d0ca4206a66be7ff6d76fc4ab8e66c6
phyloseq_to_ampvis2 <- function(physeq) {
#check object for class
if(!any(class(physeq) %in% "phyloseq"))
stop("physeq object must be of class \"phyloseq\"", call. = FALSE)
#ampvis2 requires taxonomy and abundance table, phyloseq checks for the latter
if(is.null(physeq@tax_table))
stop("No taxonomy found in the phyloseq object and is required for ampvis2", call. = FALSE)
#OTUs must be in rows, not columns
if(phyloseq::taxa_are_rows(physeq))
abund <- as.data.frame(phyloseq::otu_table(physeq)@.Data)
else
abund <- as.data.frame(t(phyloseq::otu_table(physeq)@.Data))
#tax_table is assumed to have OTUs in rows too
tax <- phyloseq::tax_table(physeq)@.Data
#merge by rownames (OTUs)
otutable <- merge(
abund,
tax,
by = 0,
all.x = TRUE,
all.y = FALSE,
sort = FALSE
)
colnames(otutable)[1] <- "OTU"
#extract sample_data (metadata)
if(!is.null(physeq@sam_data)) {
metadata <- data.frame(
phyloseq::sample_data(physeq),
row.names = phyloseq::sample_names(physeq),
stringsAsFactors = FALSE,
check.names = FALSE
)
#check if any columns match exactly with rownames
#if none matched assume row names are sample identifiers
samplesCol <- unlist(lapply(metadata, function(x) {
identical(x, rownames(metadata))}))
if(any(samplesCol)) {
#error if a column matched and it's not the first
if(!samplesCol[[1]])
stop("Sample ID's must be in the first column in the sample metadata, please reorder", call. = FALSE)
} else {
#assume rownames are sample identifiers, merge at the end with name "SampleID"
if(any(colnames(metadata) %in% "SampleID"))
stop("A column in the sample metadata is already named \"SampleID\" but does not seem to contain sample ID's", call. = FALSE)
metadata$SampleID <- rownames(metadata)
#reorder columns so SampleID is the first
metadata <- metadata[, c(which(colnames(metadata) %in% "SampleID"), 1:(ncol(metadata)-1L)), drop = FALSE]
}
} else
metadata <- NULL
#extract phylogenetic tree, assumed to be of class "phylo"
if(!is.null(physeq@phy_tree)) {
tree <- phyloseq::phy_tree(physeq)
} else
tree <- NULL
#extract OTU DNA sequences, assumed to be of class "XStringSet"
if(!is.null(physeq@refseq)) {
#convert XStringSet to DNAbin using a temporary file (easiest)
fastaTempFile <- tempfile(pattern = "ampvis2_", fileext = ".fa")
Biostrings::writeXStringSet(physeq@refseq, filepath = fastaTempFile)
} else
fastaTempFile <- NULL
#load as normally with amp_load
ampvis2::amp_load(
otutable = otutable,
metadata = metadata,
tree = tree,
fasta = fastaTempFile
)
}
#variance stabilisation from DESeq2
ps_vst <- function(ps, factor){
diagdds = phyloseq_to_deseq2(ps, as.formula(paste( "~", factor)))
diagdds = estimateSizeFactors(diagdds, type="poscounts")
diagdds = estimateDispersions(diagdds, fitType = "local")
pst <- varianceStabilizingTransformation(diagdds)
pst.dimmed <- t(as.matrix(assay(pst)))
# pst.dimmed[pst.dimmed < 0.0] <- 0.0
ps.varstab <- ps
otu_table(ps.varstab) <- otu_table(pst.dimmed, taxa_are_rows = FALSE)
return(ps.varstab)
}
#WGCNA visualisation
#result - list class object with attributes:
# ps - phyloseq object
# amp - ampwis2 object
# heat - heatmap with absalute read numbers
# heat_rel - hetmap with relative abundances
# tree - phylogenetic tree with taxonomy
color_filt <- function(ps, df){
library(tidyverse)
library(reshape2)
library(gridExtra)
l = list()
for (i in levels(df$module)){
message(i)
color_name <- df %>%
filter(module == i) %>%
pull(asv) %>%
unique()
ps.col <- prune_taxa(color_name, ps)
amp.col <- phyloseq_to_ampvis2(ps.col)
heat <- amp_heatmap(amp.col, tax_show = 60,
group_by = "Day",
tax_aggregate = "OTU",
tax_add = "Genus",
normalise=FALSE,
showRemainingTaxa = TRUE)
ps.rel <- phyloseq::transform_sample_counts(ps.col, function(x) x / sum(x) * 100)
amp.r <- phyloseq_to_ampvis2(ps.rel)
heat.rel <- amp_heatmap(amp.r, tax_show = 60,
group_by = "Day",
tax_aggregate = "OTU",
tax_add = "Genus",
normalise=FALSE,
showRemainingTaxa = TRUE)
tree <- ps.col@phy_tree
taxa <- as.data.frame(ps.col@tax_table@.Data)
p1 <- ggtree(tree) +
geom_tiplab(size=2, align=TRUE, linesize=.5) +
theme_tree2()
taxa[taxa == "Allorhizobium-Neorhizobium-Pararhizobium-Rhizobium"] <- "Allorhizobium"
taxa[taxa == "Burkholderia-Caballeronia-Paraburkholderia"] <- "Burkholderia"
tx <- taxa %>%
rownames_to_column("id") %>%
mutate(id = factor(id, levels = rev(get_taxa_name(p1)))) %>%
dplyr::select(-c(Kingdom, Species, Order)) %>%
melt(id.var = 'id')
p2 <- ggplot(tx, aes(variable, id)) +
geom_tile(aes(fill = value), alpha = 0.4) +
geom_text(aes(label = value), size = 3) +
theme_bw() +
theme(legend.position = "none",
axis.ticks.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank())
p <- ggpubr::ggarrange(p1, p2, widths = c(0.9, 1))
l[[i]] <- list("ps" = ps.col,
"amp" = amp.col,
"heat" = heat,
"heat_rel" = heat.rel,
"tree" = p,
"taxa" = knitr::kable(taxa))
}
return(l)
}
detachAllPackages <- function() {
basic.packages <- c("package:stats","package:graphics","package:grDevices","package:utils","package:datasets","package:methods","package:base")
package.list <- search()[ifelse(unlist(gregexpr("package:",search()))==1,TRUE,FALSE)]
package.list <- setdiff(package.list,basic.packages)
if (length(package.list)>0) for (package in package.list) detach(package, character.only=TRUE, force = TRUE)
}
plot_alpha_w_toc <- function(ps, group, metric) {
require(phyloseq)
require(ggplot2)
ps_a <- prune_taxa(taxa_sums(ps) > 0, ps)
er <- estimate_richness(ps_a)
df_er <- cbind(ps_a@sam_data, er)
df_er <- df_er %>% select(c(group, metric))
stat.test <- aov(as.formula(paste0(metric, "~", group)), data = df_er) %>%
rstatix::tukey_hsd()
y <- seq(max(er[[metric]]), length=length(stat.test$p.adj), by=max(er[[metric]]/20))
plot_richness(ps_a, x=group, measures=metric) +
geom_boxplot() +
geom_point(size=1.2, alpha=0.3) +
ggpubr::stat_pvalue_manual(
stat.test,
label = "p.adj.signif",
y.position = y) +
theme_light() +
scale_color_brewer(palette="Dark2") +
theme(axis.text.x = element_text(angle = 90),
axis.title.x=element_blank()) +
labs(y=paste(metric, "index"))
}
# standart NMDS plot tool frop phyloseq with some additional aestatics
# have stress value on plot - may work as fuck
beta_custom_norm_NMDS_elli_w <- function(ps, seed = 7888, normtype="vst", Color="What", Group="Repeat"){
require(phyloseq)
require(ggplot2)
require(ggpubr)
library(ggforce)
ordination.b <- ordinate(ps, "NMDS", "bray")
mds <- as.data.frame(ordination.b$points)
p <- plot_ordination(ps,
ordination.b,
type="sample",
color = Color,
title="NMDS - Bray-Curtis",
# title=NULL,
axes = c(2,3) ) +
theme_bw() +
theme(text = element_text(size = 10)) +
geom_point(size = 3) +
annotate("text",
x=min(mds$MDS1) + abs(min(mds$MDS1))/7,
y=max(mds$MDS2),
label=paste0("Stress -- ", round(ordination.b$stress, 3))) +
geom_mark_ellipse(aes_string(group = Group, label = Group),
label.fontsize = 10,
label.buffer = unit(2, "mm"),
label.minwidth = unit(5, "mm"),
con.cap = unit(0.1, "mm"),
con.colour='gray') +
theme(legend.position = "none") +
ggplot2::scale_fill_viridis_c(option = "H")
return(p)
}
# alpha with aov + tukie post-hock - useless, but it looks pretty good
plot_alpha_w_toc <- function(ps, group, metric) {
require(phyloseq)
require(ggplot2)
ps_a <- prune_taxa(taxa_sums(ps) > 0, ps)
er <- estimate_richness(ps_a)
df_er <- cbind(ps_a@sam_data, er)
df_er <- df_er %>% select(c(group, metric))
stat.test <- aov(as.formula(paste0(metric, "~", group)), data = df_er) %>%
rstatix::tukey_hsd()
y <- seq(max(er[[metric]]), length=length(stat.test$p.adj.signif[stat.test$p.adj.signif != "ns"]), by=max(er[[metric]]/20))
plot_richness(ps_a, x=group, measures=metric) +
geom_boxplot() +
geom_point(size=1.2, alpha=0.3) +
stat_pvalue_manual(
stat.test,
label = "p.adj.signif",
y.position = y,
hide.ns=TRUE) +
theme_light() +
scale_color_brewer(palette="Dark2") +
theme(axis.text.x = element_text(angle = 90),
axis.title.x=element_blank()) +
labs(y=paste(metric, "index"))
}
p1 <- plot_alpha_w_toc(ps = ps.f, group = "Day", metric = "Observed")
## Note: Using an external vector in selections is ambiguous.
## ℹ Use `all_of(group)` instead of `group` to silence this message.
## ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This message is displayed once per session.
## Note: Using an external vector in selections is ambiguous.
## ℹ Use `all_of(metric)` instead of `metric` to silence this message.
## ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This message is displayed once per session.
p2 <- plot_alpha_w_toc(ps = ps.f, group = "Day", metric = "Shannon")
p3 <- plot_alpha_w_toc(ps = ps.f, group = "Day", metric = "InvSimpson")
ggarrange(p1, p2, p3, nrow = 1)
Add Group parameter to metadata -
- early - D01, D03, D05
- middle - D07, D08, D10
- late - D13, D14, D15
sample.data <- ps.f@sam_data %>%
data.frame() %>%
mutate(Group = if_else(Day %in% c("D01", "D03", "D05"), "early",
if_else(Day %in% c("D07", "D08","D10"), "middle", "late"))) %>%
mutate(Group = factor(Group, levels=c("early", "middle","late"))) %>%
phyloseq::sample_data()
sample_data(ps.f) <- sample.data
p.observed <- plot_alpha_w_toc(ps = ps.f, group = "Group", metric = c("Observed")) +
theme(axis.title.y = element_blank())
p.shannon <- plot_alpha_w_toc(ps = ps.f, group = "Group", metric = c("Shannon")) +
theme(axis.title.y = element_blank())
p.simpson <- plot_alpha_w_toc(ps = ps.f, group = "Group", metric = c("InvSimpson")) +
theme(axis.title.y = element_blank())
ggpubr::ggarrange(p.observed, p.shannon, p.simpson, ncol = 3)
permanova - group significantlly different - dispersion between is more, than inside groups
dist <- phyloseq::distance(ps.f, "bray")
metadata <- as(sample_data(ps.f@sam_data), "data.frame")
vegan::adonis2(dist ~ Group, data = metadata)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## vegan::adonis2(formula = dist ~ Group, data = metadata)
## Df SumOfSqs R2 F Pr(>F)
## Group 2 3.2742 0.33893 8.2033 0.001 ***
## Residual 32 6.3861 0.66107
## Total 34 9.6604 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Дистанция брей кертиса по шагам между каждыми днями. (Например - все D15 против D14) Может на эту картинку подвесить еще почвенное дыхание?
ps.f.r <- rarefy_even_depth(ps.f, rngseed = 777)
avg.r <- ps.f.r@otu_table %>%
as.data.frame() %>%
vegan::avgdist(10)
avg <- ps.f@otu_table %>%
as.data.frame() %>%
vegan::avgdist(10)
# avg %>%
# as.matrix() %>%
# as_tibble(rownames= "sample") %>%
# pivot_longer(-sample) %>%
# filter(sample < name) %>%
# mutate(repeat_a = str_replace(sample, ".*-", ""),
# repeat_b = str_replace(name, ".*-", ""),
# day_a = as.numeric(str_replace(sapply(strsplit(sample, "-"), `[`, 3), "D", "")),
# day_b = as.numeric(str_replace(sapply(strsplit(name, "-"), `[`, 3), "D", "")),
# diff = abs(day_a - day_b),
# early = day_a < 10) %>%
# filter(repeat_a == repeat_b & diff < 10) %>%
# group_by(diff, repeat_a, early) %>%
# summarize(median = median(value)) %>%
# ungroup() %>%
# ggplot(aes(x=diff, y=median, color=early, group=paste0(repeat_a, early))) +
# geom_line(size=0.25) +
# geom_smooth(aes(group=early), se=FALSE, size=4) +
# labs(x="Distance between time points",
# y="Median Bray-Curtis distance") +
# scale_x_continuous(breaks=1:9) +
# scale_color_manual(name=NULL,
# breaks=c(TRUE, FALSE),
# values=c("blue", "red"),
# labels=c("Early", "Late")) +
# guides(color = guide_legend(override.aes = list(size=1))) +
# theme_classic()
avg.r %>%
as.matrix() %>%
as_tibble(rownames= "sample") %>%
pivot_longer(-sample) %>%
filter(sample < name) %>%
mutate(repeat_a = str_replace(sample, ".*-", ""),
repeat_b = str_replace(name, ".*-", ""),
day_a = as.numeric(str_replace(sapply(strsplit(sample, "-"), `[`, 3), "D", "")),
day_b = as.numeric(str_replace(sapply(strsplit(name, "-"), `[`, 3), "D", ""))) %>%
mutate(day_a = as.numeric(as.factor(day_a) %>% forcats::fct_recode("2" = "3", "3" = "5", "4" = "7", "5" = "8", "6" = "10", "7" = "13", "8" = "14", "9" = "14", "10" = "15")),
day_b = as.numeric(as.factor(day_b) %>% forcats::fct_recode("2" = "3", "3" = "5", "4" = "7", "5" = "8", "6" = "10", "7" = "12", "8" = "13", "9" = "14", "10" = "15")),
diff = abs(day_a - day_b)) %>%
filter(diff == 1) %>%
mutate(day_b = as.factor(day_b) %>% forcats::fct_recode("D03" = "2", "D05" = "3","D07" = "4", "D08" = "5", "D10" = "6", "D13" = "7", "D14" = "8", "D15" = "10", "D12" = "7")) %>%
ggplot(aes(x=day_b, y=value)) +
geom_boxplot() +
theme_bw() +
labs(x="Time points",
y="Bray-Curtis distance")
IN WORK!!
beta_custom_norm_NMDS_elli_w(ps.f, C="Group", G="Day")
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1048129
## Run 1 stress 0.1372731
## Run 2 stress 0.1318495
## Run 3 stress 0.1048129
## ... Procrustes: rmse 0.00004856956 max resid 0.0001949092
## ... Similar to previous best
## Run 4 stress 0.1048129
## ... New best solution
## ... Procrustes: rmse 0.00001842645 max resid 0.00006948497
## ... Similar to previous best
## Run 5 stress 0.1426236
## Run 6 stress 0.1048204
## ... Procrustes: rmse 0.001770292 max resid 0.007527799
## ... Similar to previous best
## Run 7 stress 0.1064532
## Run 8 stress 0.1064733
## Run 9 stress 0.1048204
## ... Procrustes: rmse 0.001767644 max resid 0.007510715
## ... Similar to previous best
## Run 10 stress 0.1064532
## Run 11 stress 0.1318495
## Run 12 stress 0.1064733
## Run 13 stress 0.1064733
## Run 14 stress 0.1460032
## Run 15 stress 0.1048204
## ... Procrustes: rmse 0.001769848 max resid 0.007525048
## ... Similar to previous best
## Run 16 stress 0.1064532
## Run 17 stress 0.1372728
## Run 18 stress 0.1048129
## ... Procrustes: rmse 0.000002711789 max resid 0.0000100797
## ... Similar to previous best
## Run 19 stress 0.1460029
## Run 20 stress 0.146003
## *** Solution reached
Разделим датасет на две группы
1-я - в более чем 10% образцов должно быть хотя бы 10 ридов - эта группа пойдет в анализ далее
ps.inall <- phyloseq::filter_taxa(ps.f, function(x) sum(x > 10) > (0.1*length(x)), TRUE)
amp.inall <- phyloseq_to_ampvis2(ps.inall)
amp_heatmap(amp.inall,
tax_show = 40,
group_by = "Day",
tax_aggregate = "OTU",
tax_add = "Genus",
normalise=FALSE,
showRemainingTaxa = TRUE)
Вторая - оставшиеся, но более 100 ридов по всем образцам(далее -
вылетающие мажоры(ВМ))
Остальные филотипы выкидываем из анализа
ps.exl <- phyloseq::filter_taxa(ps.f, function(x) sum(x > 10) < (0.1*length(x)), TRUE)
ps.exl <- prune_taxa(taxa_sums(ps.exl) > 100, ps.exl)
amp.exl <- phyloseq_to_ampvis2(ps.exl)
amp_heatmap(amp.exl,
tax_show = 40,
group_by = "Day",
tax_aggregate = "OTU",
tax_add = "Genus",
normalise=FALSE,
showRemainingTaxa = TRUE)
То же, но c относительной представленностью.
ps.per <- phyloseq::transform_sample_counts(ps.f, function(x) x / sum(x) * 100)
ps.exl.taxa <- taxa_names(ps.exl)
ps.per.exl <- prune_taxa(ps.exl.taxa, ps.per)
amp.exl.r <- phyloseq_to_ampvis2(ps.per.exl)
amp_heatmap(amp.exl.r, tax_show = 60,
group_by = "Day",
tax_aggregate = "OTU",
tax_add = "Genus",
round = 2,
normalise=FALSE,
showRemainingTaxa = TRUE)
amp_heatmap(amp.exl.r, tax_show = 60,
group_by = "Day",
tax_aggregate = "OTU",
tax_add = "Genus",
round = 2,
normalise=FALSE, )
ВМ объединенные по родам. Относительная представленность.
amp_heatmap(amp.exl.r,
tax_show = 30,
group_by = "Day",
tax_aggregate = "Genus",
tax_add = "Phylum",
tax_class = "Proteobacteria",
round = 2,
normalise=FALSE,
showRemainingTaxa = TRUE)
Если посмотреть, если как распределяются ВМ по дням - похоже распределение ВМ связано не только с особенностью отдельных мешочков, но и с чисто техническими особенностями - вылетающие значения вылетают в основном не с биологическими повторностями, а с техническими(красная линия - выбрана визуально)
p_box <- phyloseq::sample_sums(ps.per.exl) %>%
as.data.frame(col.names = "values") %>%
setNames(., nm = "values") %>%
rownames_to_column("samples") %>%
mutate(Day = sapply(strsplit(samples, "-"), `[`, 3)) %>%
ggplot(aes(x=Day, y=values, color=Day, fill = Day)) +
geom_boxplot(aes(color=Day, fill = Day)) +
geom_point(color = "black", position = position_dodge(width=0.2)) +
geom_hline(yintercept = 10, colour = "red") +
theme_bw() +
theme(legend.position = "none")
p_box <- p_box + viridis::scale_color_viridis(option = "H", discrete = TRUE, direction=1, begin=0.1, end = 0.9, alpha = 0.5)
p_box + viridis::scale_fill_viridis(option = "H", discrete = TRUE, direction=1, begin=0.1, end = 0.9, alpha = 0.3)
после clr нормализации / выглядит отвратительно попробуем нормализацию vst из DESeq2
otu.inall <- as.data.frame(ps.inall@otu_table@.Data)
ps.inall.clr <- ps.inall
otu_table(ps.inall.clr) <- phyloseq::otu_table(compositions::clr(otu.inall), taxa_are_rows = FALSE)
data <- ps.inall.clr@otu_table@.Data %>%
as.data.frame()
rownames(data) <- as.character(ps.inall.clr@sam_data$Description)
powers <- c(c(1:10), seq(from = 12, to=30, by=1))
sft <- pickSoftThreshold(data, powerVector = powers, verbose = 5, networkType = "signed")
## pickSoftThreshold: will use block size 338.
## pickSoftThreshold: calculating connectivity for given powers...
## ..working on genes 1 through 338 of 338
## Warning: executing %dopar% sequentially: no parallel backend registered
## Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k.
## 1 1 0.185 23.60 0.8790 168.0000 168.00000 172.000
## 2 2 0.199 -15.50 0.8970 87.6000 87.20000 95.100
## 3 3 0.518 -13.10 0.8580 47.3000 46.80000 55.900
## 4 4 0.390 -6.46 0.8050 26.5000 26.00000 34.500
## 5 5 0.262 -3.46 0.7820 15.4000 14.90000 22.300
## 6 6 0.232 -2.30 0.7770 9.2300 8.78000 15.000
## 7 7 0.251 -1.87 0.6990 5.7300 5.36000 10.500
## 8 8 0.319 -1.75 0.7570 3.6800 3.41000 7.620
## 9 9 0.527 -2.07 0.7770 2.4400 2.20000 5.890
## 10 10 0.648 -2.28 0.8490 1.6600 1.47000 4.690
## 11 12 0.827 -2.41 0.8970 0.8420 0.69300 3.190
## 12 13 0.793 -2.38 0.7930 0.6220 0.49000 2.700
## 13 14 0.277 -3.97 0.0914 0.4690 0.35500 2.320
## 14 15 0.282 -3.83 0.1000 0.3610 0.26100 2.010
## 15 16 0.894 -2.17 0.9030 0.2820 0.19800 1.770
## 16 17 0.905 -2.10 0.9080 0.2240 0.15000 1.560
## 17 18 0.922 -2.04 0.9240 0.1810 0.11500 1.390
## 18 19 0.927 -2.01 0.9230 0.1480 0.08820 1.250
## 19 20 0.938 -1.93 0.9370 0.1220 0.06780 1.130
## 20 21 0.946 -1.86 0.9470 0.1020 0.05260 1.020
## 21 22 0.315 -3.12 0.2220 0.0856 0.04140 0.926
## 22 23 0.303 -3.96 0.1560 0.0728 0.03270 0.844
## 23 24 0.307 -2.92 0.1940 0.0623 0.02590 0.772
## 24 25 0.303 -2.84 0.1800 0.0538 0.02060 0.708
## 25 26 0.921 -1.63 0.8990 0.0467 0.01650 0.657
## 26 27 0.892 -1.61 0.8650 0.0408 0.01320 0.622
## 27 28 0.155 -1.93 -0.0202 0.0359 0.01040 0.591
## 28 29 0.199 -2.77 -0.0228 0.0317 0.00841 0.563
## 29 30 0.199 -2.69 -0.0203 0.0282 0.00684 0.537
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2], xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n", main = paste("Scale independence"))
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2], labels=powers,cex=0.9,col="red")
abline(h=0.9,col="salmon")
после vst нормализации
ps.varstab <- ps_vst(ps.inall, "Day")
data2 <- ps.varstab@otu_table@.Data %>%
as.data.frame()
rownames(data2) <- as.character(ps.varstab@sam_data$Description)
powers <- c(seq(from = 1, to=10, by=0.5), seq(from = 11, to=20, by=1))
sft2 <- pickSoftThreshold(data2, powerVector = powers, verbose = 5, networkType = "signed hybrid")
## pickSoftThreshold: will use block size 338.
## pickSoftThreshold: calculating connectivity for given powers...
## ..working on genes 1 through 338 of 338
## Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k.
## 1 1.0 0.409 -1.05 0.704 50.100 45.60000 97.80
## 2 1.5 0.578 -1.07 0.895 30.900 27.00000 71.20
## 3 2.0 0.723 -1.15 0.942 20.300 17.90000 53.60
## 4 2.5 0.767 -1.24 0.871 14.000 12.00000 41.70
## 5 3.0 0.809 -1.31 0.909 10.000 7.86000 33.40
## 6 3.5 0.872 -1.31 0.958 7.380 5.39000 27.30
## 7 4.0 0.888 -1.30 0.968 5.580 3.76000 22.60
## 8 4.5 0.858 -1.36 0.910 4.320 2.67000 19.00
## 9 5.0 0.892 -1.36 0.943 3.400 1.94000 16.20
## 10 5.5 0.899 -1.38 0.954 2.720 1.41000 13.90
## 11 6.0 0.896 -1.39 0.962 2.210 1.07000 12.10
## 12 6.5 0.897 -1.33 0.962 1.820 0.84700 10.60
## 13 7.0 0.910 -1.30 0.970 1.520 0.65400 9.30
## 14 7.5 0.914 -1.31 0.966 1.280 0.51500 8.25
## 15 8.0 0.913 -1.28 0.954 1.090 0.41500 7.36
## 16 8.5 0.905 -1.28 0.937 0.938 0.34000 6.59
## 17 9.0 0.895 -1.30 0.913 0.813 0.27800 5.99
## 18 9.5 0.878 -1.33 0.900 0.709 0.23200 5.49
## 19 10.0 0.886 -1.34 0.894 0.624 0.19100 5.07
## 20 11.0 0.911 -1.31 0.928 0.491 0.12500 4.35
## 21 12.0 0.890 -1.29 0.894 0.396 0.08590 3.79
## 22 13.0 0.910 -1.26 0.918 0.325 0.06120 3.33
## 23 14.0 0.882 -1.25 0.856 0.272 0.04420 2.95
## 24 15.0 0.889 -1.19 0.864 0.230 0.03230 2.63
## 25 16.0 0.913 -1.18 0.893 0.198 0.02240 2.35
## 26 17.0 0.930 -1.16 0.912 0.172 0.01620 2.20
## 27 18.0 0.955 -1.19 0.947 0.151 0.01190 2.09
## 28 19.0 0.908 -1.21 0.896 0.134 0.00896 1.99
## 29 20.0 0.916 -1.21 0.906 0.120 0.00649 1.91
plot(sft2$fitIndices[,1], -sign(sft2$fitIndices[,3])*sft2$fitIndices[,2], xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n", main = paste("Scale independence"))
text(sft2$fitIndices[,1], -sign(sft2$fitIndices[,3])*sft2$fitIndices[,2], labels=powers,cex=0.9,col="red")
abline(h=0.9,col="salmon")
# WGCNA, as old and not supported
detachAllPackages()
library(WGCNA)
net3 <- WGCNA::blockwiseModules(data2,
power=5.5,
TOMType="signed",
networkType="signed hybrid",
nThreads=0)
mergedColors2 <- WGCNA::labels2colors(net3$colors, colorSeq = c("salmon", "darkgreen", "cyan", "red", "blue", "plum"))
plotDendroAndColors(
net3$dendrograms[[1]],
mergedColors2[net3$blockGenes[[1]]],
"Module colors",
dendroLabels = FALSE,
hang = 0.03,
addGuide = TRUE,
guideHang = 0.05)
library(phyloseq)
library(tidyverse)
library(ggpubr)
library(ampvis2)
library(heatmaply)
library(WGCNA)
library(phyloseq)
library(ggtree)
library(tidyverse)
library(KneeArrower)
modules_of_interest = mergedColors2 %>%
unique()
module_df <- data.frame(
asv = names(net3$colors),
colors = mergedColors2
)
# module_df[module_df == "yellow"] <- "salmon"
submod <- module_df %>%
subset(colors %in% modules_of_interest)
row.names(module_df) = module_df$asv
subexpr = as.data.frame(t(data2))[submod$asv,]
submod_df <- data.frame(subexpr) %>%
mutate(
asv = row.names(.)
) %>%
pivot_longer(-asv) %>%
mutate(
module = module_df[asv,]$colors
)
submod_df <- submod_df %>%
mutate(name = gsub("\\_.*","",submod_df$name)) %>%
group_by(name, asv) %>%
summarise(value = mean(value), asv = asv, module = module) %>%
relocate(c(asv, name, value, module)) %>%
ungroup() %>%
mutate(module = as.factor(module))
p <- submod_df %>%
ggplot(., aes(x=name, y=value, group=asv)) +
geom_line(aes(color = module),
alpha = 0.2) +
theme_bw() +
theme(
axis.text.x = element_text(angle = 90),
legend.position = "none") +
facet_grid(rows = vars(module)) +
labs(x = "treatment",
y = "normalized expression")
p + scale_color_manual(values = levels(submod_df$module))
UNIFRAC - колбаска сужается
ps.inall.col <- ps.inall
df <- module_df %>%
rename("id" = "asv")
df <- df %>%
dplyr::select(-"id") %>%
mutate(colors = as.factor(colors))
taxa <- as.data.frame(ps.inall@tax_table@.Data)
tx <- cbind(taxa, df)
tx$colors <- factor(tx$colors, levels = c("salmon", "darkgreen", "cyan", "red", "blue", "plum"))
tax_table(ps.inall.col) <- tax_table(as.matrix(tx))
ord <- ordinate(ps.inall.col, "NMDS", "unifrac")
## Run 0 stress 0.0685455
## Run 1 stress 0.05568715
## ... New best solution
## ... Procrustes: rmse 0.04076806 max resid 0.2160826
## Run 2 stress 0.06784429
## Run 3 stress 0.05628488
## Run 4 stress 0.05568719
## ... Procrustes: rmse 0.0001769095 max resid 0.0006062451
## ... Similar to previous best
## Run 5 stress 0.07937976
## Run 6 stress 0.07544596
## Run 7 stress 0.05628489
## Run 8 stress 0.05568719
## ... Procrustes: rmse 0.00002559041 max resid 0.00008856774
## ... Similar to previous best
## Run 9 stress 0.05568712
## ... New best solution
## ... Procrustes: rmse 0.00002106782 max resid 0.00006934972
## ... Similar to previous best
## Run 10 stress 0.06740137
## Run 11 stress 0.06740136
## Run 12 stress 0.06784421
## Run 13 stress 0.06784436
## Run 14 stress 0.06854546
## Run 15 stress 0.07953205
## Run 16 stress 0.05568719
## ... Procrustes: rmse 0.00004255519 max resid 0.0001401031
## ... Similar to previous best
## Run 17 stress 0.06740138
## Run 18 stress 0.06728216
## Run 19 stress 0.05568718
## ... Procrustes: rmse 0.0000373923 max resid 0.0001312747
## ... Similar to previous best
## Run 20 stress 0.07888487
## *** Solution reached
plot_ordination(ps.inall.col, ord, type = "species", color = "colors") +
scale_color_manual(values = c("salmon", "darkgreen", "cyan", "red", "blue", "plum")) +
theme_bw() +
theme(legend.position = "none")
Bray - колбаска равномерна
Поздние стадии слева - при этом ранние кластеры накладываются, поздние
разделены Что влияет на ось2? Явно есть какой-то паттерн.
ord <- ordinate(ps.inall.col, "NMDS", "bray")
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.09036937
## Run 1 stress 0.0903694
## ... Procrustes: rmse 0.00001927814 max resid 0.00007086537
## ... Similar to previous best
## Run 2 stress 0.09024752
## ... New best solution
## ... Procrustes: rmse 0.003386048 max resid 0.01505777
## Run 3 stress 0.09036938
## ... Procrustes: rmse 0.003390695 max resid 0.01505826
## Run 4 stress 0.09036937
## ... Procrustes: rmse 0.003386035 max resid 0.01507802
## Run 5 stress 0.09036937
## ... Procrustes: rmse 0.003386432 max resid 0.01505862
## Run 6 stress 0.1257473
## Run 7 stress 0.09036937
## ... Procrustes: rmse 0.003385691 max resid 0.01507976
## Run 8 stress 0.1268919
## Run 9 stress 0.09036944
## ... Procrustes: rmse 0.003410952 max resid 0.01523599
## Run 10 stress 0.09024752
## ... New best solution
## ... Procrustes: rmse 0.00003844186 max resid 0.0001638862
## ... Similar to previous best
## Run 11 stress 0.09036937
## ... Procrustes: rmse 0.003380784 max resid 0.0150593
## Run 12 stress 0.09036938
## ... Procrustes: rmse 0.003379871 max resid 0.01508219
## Run 13 stress 0.09036939
## ... Procrustes: rmse 0.003385498 max resid 0.01512324
## Run 14 stress 0.1358024
## Run 15 stress 0.09024754
## ... Procrustes: rmse 0.00003794727 max resid 0.0001665169
## ... Similar to previous best
## Run 16 stress 0.1258672
## Run 17 stress 0.1358024
## Run 18 stress 0.09036937
## ... Procrustes: rmse 0.003376715 max resid 0.01503468
## Run 19 stress 0.126985
## Run 20 stress 0.1360016
## *** Solution reached
plot_ordination(ps.inall.col, ord, type = "species", color = "colors") +
scale_color_manual(values = c("salmon", "darkgreen", "cyan", "red", "blue", "plum")) +
theme_bw() +
theme(legend.position = "none")
Далее идут одинаковые картинки по всем группам.
l_vst <- color_filt(ps.inall, submod_df)
l_vst
$blue \(blue\)ps phyloseq-class experiment-level object otu_table() OTU Table: [ 86 taxa and 35 samples ] sample_data() Sample Data: [ 35 samples by 3 sample variables ] tax_table() Taxonomy Table: [ 86 taxa by 7 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 86 tips and 85 internal nodes ] refseq() DNAStringSet: [ 86 reference sequences ]
\(blue\)amp ampvis2 object with 5 elements. Summary of OTU table: Samples OTUs Total#Reads Min#Reads Max#Reads Median#Reads 35 86 89493 0 11387 1819 Avg#Reads 2556.94
Assigned taxonomy: Kingdom Phylum Class Order Family Genus Species 86(100%) 86(100%) 86(100%) 85(98.84%) 79(91.86%) 50(58.14%) 7(8.14%)
Metadata variables: 4 SampleID, Day, Description, Group
\(blue\)heat
\(blue\)heat_rel
\(blue\)tree
\(blue\)taxa
| Kingdom | Phylum | Class | Order | Family | Genus | Species | |
|---|---|---|---|---|---|---|---|
| Seq314 | Bacteria | Cyanobacteria | Vampirivibrionia | Obscuribacterales | Obscuribacteraceae | NA | NA |
| Seq23 | Bacteria | Cyanobacteria | Vampirivibrionia | Obscuribacterales | Obscuribacteraceae | NA | NA |
| Seq65 | Bacteria | Proteobacteria | Alphaproteobacteria | Sphingomonadales | Sphingomonadaceae | Sphingomonas | NA |
| Seq188 | Bacteria | Proteobacteria | Alphaproteobacteria | Reyranellales | Reyranellaceae | Reyranella | soli |
| Seq334 | Bacteria | Proteobacteria | Alphaproteobacteria | Reyranellales | Reyranellaceae | Reyranella | NA |
| Seq75 | Bacteria | Proteobacteria | Alphaproteobacteria | Reyranellales | Reyranellaceae | Reyranella | NA |
| Seq167 | Bacteria | Proteobacteria | Alphaproteobacteria | Reyranellales | Reyranellaceae | Reyranella | massiliensis |
| Seq273 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhizobiales | Kaistiaceae | Kaistia | NA |
| Seq131 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhizobiales | Xanthobacteraceae | Pseudolabrys | NA |
| Seq138 | Bacteria | Proteobacteria | Alphaproteobacteria | Micropepsales | Micropepsaceae | NA | NA |
| Seq104 | Bacteria | Proteobacteria | Alphaproteobacteria | Micropepsales | Micropepsaceae | NA | NA |
| Seq129 | Bacteria | Proteobacteria | Alphaproteobacteria | Sphingomonadales | Sphingomonadaceae | Altererythrobacter | NA |
| Seq28 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhizobiales | Devosiaceae | Devosia | NA |
| Seq360 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhizobiales | Rhizobiaceae | Mesorhizobium | NA |
| Seq292 | Bacteria | Proteobacteria | Alphaproteobacteria | Caulobacterales | Caulobacteraceae | NA | NA |
| Seq166 | Bacteria | Proteobacteria | Alphaproteobacteria | Caulobacterales | Hyphomonadaceae | Hirschia | NA |
| Seq39 | Bacteria | Gemmatimonadota | Gemmatimonadetes | Gemmatimonadales | Gemmatimonadaceae | NA | NA |
| Seq30 | Bacteria | Actinobacteriota | Thermoleophilia | Solirubrobacterales | Solirubrobacteraceae | Conexibacter | NA |
| Seq151 | Bacteria | Actinobacteriota | Thermoleophilia | Solirubrobacterales | 67-14 | NA | NA |
| Seq250 | Bacteria | Actinobacteriota | Thermoleophilia | Solirubrobacterales | 67-14 | NA | NA |
| Seq135 | Bacteria | Myxococcota | Polyangia | Polyangiales | BIrii41 | NA | NA |
| Seq46 | Bacteria | Myxococcota | Polyangia | Polyangiales | BIrii41 | NA | NA |
| Seq35 | Bacteria | Myxococcota | Polyangia | Polyangiales | BIrii41 | NA | NA |
| Seq106 | Bacteria | Myxococcota | Polyangia | Polyangiales | Sandaracinaceae | NA | NA |
| Seq345 | Bacteria | Bdellovibrionota | Oligoflexia | 0319-6G20 | NA | NA | NA |
| Seq198 | Bacteria | Bdellovibrionota | Bdellovibrionia | Bdellovibrionales | Bdellovibrionaceae | Bdellovibrio | NA |
| Seq97 | Bacteria | Bdellovibrionota | Bdellovibrionia | Bdellovibrionales | Bdellovibrionaceae | Bdellovibrio | NA |
| Seq183 | Bacteria | Dependentiae | Babeliae | Babeliales | UBA12409 | NA | NA |
| Seq424 | Bacteria | Actinobacteriota | Actinobacteria | Propionibacteriales | Nocardioidaceae | Kribbella | NA |
| Seq136 | Bacteria | Actinobacteriota | Acidimicrobiia | Microtrichales | Iamiaceae | Iamia | NA |
| Seq174 | Bacteria | Actinobacteriota | Acidimicrobiia | Microtrichales | Iamiaceae | Iamia | NA |
| Seq43 | Bacteria | Actinobacteriota | Actinobacteria | Micrococcales | Microbacteriaceae | Galbitalea | NA |
| Seq420 | Bacteria | Actinobacteriota | Actinobacteria | Micromonosporales | Micromonosporaceae | Luedemannella | NA |
| Seq127 | Bacteria | Actinobacteriota | Actinobacteria | Micromonosporales | Micromonosporaceae | Dactylosporangium | NA |
| Seq288 | Bacteria | Actinobacteriota | Actinobacteria | Streptomycetales | Streptomycetaceae | NA | NA |
| Seq449 | Bacteria | Actinobacteriota | Actinobacteria | Propionibacteriales | Propionibacteriaceae | Jiangella | NA |
| Seq91 | Bacteria | Bacteroidota | Bacteroidia | Bacteroidales | NA | NA | NA |
| Seq149 | Bacteria | Bacteroidota | Bacteroidia | Bacteroidetes VC2.1 Bac22 | NA | NA | NA |
| Seq169 | Bacteria | Bacteroidota | Bacteroidia | Sphingobacteriales | Sphingobacteriaceae | Mucilaginibacter | calamicampi |
| Seq81 | Bacteria | Bacteroidota | Bacteroidia | Sphingobacteriales | NS11-12 marine group | NA | NA |
| Seq211 | Bacteria | Bacteroidota | Bacteroidia | NA | NA | NA | NA |
| Seq119 | Bacteria | Bacteroidota | Bacteroidia | Sphingobacteriales | env.OPS 17 | NA | NA |
| Seq77 | Bacteria | Spirochaetota | Spirochaetia | Spirochaetales | Spirochaetaceae | Salinispira | NA |
| Seq175 | Bacteria | Spirochaetota | Spirochaetia | Spirochaetales | Spirochaetaceae | Spirochaeta 2 | NA |
| Seq305 | Bacteria | Chloroflexi | Anaerolineae | Anaerolineales | Anaerolineaceae | NA | NA |
| Seq379 | Bacteria | Chloroflexi | Chloroflexia | Chloroflexales | Roseiflexaceae | NA | NA |
| Seq165 | Bacteria | Chloroflexi | Chloroflexia | Chloroflexales | Roseiflexaceae | NA | NA |
| Seq327 | Bacteria | Chloroflexi | Chloroflexia | Thermomicrobiales | JG30-KF-CM45 | NA | NA |
| Seq252 | Bacteria | Armatimonadota | Fimbriimonadia | Fimbriimonadales | Fimbriimonadaceae | NA | NA |
| Seq382 | Bacteria | Verrucomicrobiota | Verrucomicrobiae | Pedosphaerales | Pedosphaeraceae | NA | NA |
| Seq123 | Bacteria | Verrucomicrobiota | Verrucomicrobiae | Opitutales | Opitutaceae | Lacunisphaera | limnophila |
| Seq124 | Bacteria | Verrucomicrobiota | Verrucomicrobiae | Chthoniobacterales | Terrimicrobiaceae | Terrimicrobium | NA |
| Seq325 | Bacteria | Verrucomicrobiota | Verrucomicrobiae | Chthoniobacterales | Chthoniobacteraceae | LD29 | NA |
| Seq272 | Bacteria | Verrucomicrobiota | Verrucomicrobiae | Verrucomicrobiales | Verrucomicrobiaceae | Roseimicrobium | gellanilyticum |
| Seq402 | Bacteria | Verrucomicrobiota | Verrucomicrobiae | Verrucomicrobiales | Rubritaleaceae | Luteolibacter | NA |
| Seq322 | Bacteria | Planctomycetota | Planctomycetes | Pirellulales | Pirellulaceae | NA | NA |
| Seq356 | Bacteria | Planctomycetota | Planctomycetes | Pirellulales | Pirellulaceae | NA | NA |
| Seq259 | Bacteria | Planctomycetota | Planctomycetes | Pirellulales | Pirellulaceae | Pir4 lineage | NA |
| Seq279 | Bacteria | Planctomycetota | Planctomycetes | Pirellulales | Pirellulaceae | Pir4 lineage | NA |
| Seq214 | Bacteria | Planctomycetota | Planctomycetes | Planctomycetales | Rubinisphaeraceae | SH-PL14 | NA |
| Seq209 | Bacteria | Planctomycetota | Planctomycetes | Planctomycetales | Schlesneriaceae | Schlesneria | NA |
| Seq229 | Bacteria | Acidobacteriota | Vicinamibacteria | Vicinamibacterales | Vicinamibacteraceae | NA | NA |
| Seq275 | Bacteria | Acidobacteriota | Vicinamibacteria | Vicinamibacterales | Vicinamibacteraceae | NA | NA |
| Seq247 | Bacteria | Acidobacteriota | Vicinamibacteria | Vicinamibacterales | NA | NA | NA |
| Seq316 | Bacteria | Bacteroidota | Kapabacteria | Kapabacteriales | NA | NA | NA |
| Seq371 | Bacteria | Bacteroidota | Bacteroidia | Cytophagales | Microscillaceae | Ohtaekwangia | NA |
| Seq12 | Bacteria | Bacteroidota | Bacteroidia | Cytophagales | Microscillaceae | Ohtaekwangia | NA |
| Seq153 | Bacteria | Bacteroidota | Bacteroidia | Cytophagales | Microscillaceae | Ohtaekwangia | NA |
| Seq326 | Bacteria | Bacteroidota | Bacteroidia | Cytophagales | Microscillaceae | Ohtaekwangia | NA |
| Seq10 | Bacteria | Bacteroidota | Bacteroidia | Cytophagales | Microscillaceae | NA | NA |
| Seq82 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Terrimonas | NA |
| Seq462 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Pseudoflavitalea | NA |
| Seq150 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Flavitalea | NA |
| Seq290 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Edaphobaculum | NA |
| Seq339 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Edaphobaculum | NA |
| Seq120 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Taibaiella | NA |
| Seq233 | Bacteria | Proteobacteria | Gammaproteobacteria | Pseudomonadales | 211ds20 | NA | NA |
| Seq133 | Bacteria | Proteobacteria | Gammaproteobacteria | Steroidobacterales | Steroidobacteraceae | Steroidobacter | flavus |
| Seq349 | Bacteria | Proteobacteria | Gammaproteobacteria | Steroidobacterales | Steroidobacteraceae | Steroidobacter | NA |
| Seq315 | Bacteria | Proteobacteria | Gammaproteobacteria | Gammaproteobacteria Incertae Sedis | Unknown Family | Acidibacter | NA |
| Seq57 | Bacteria | Proteobacteria | Gammaproteobacteria | Gammaproteobacteria Incertae Sedis | Unknown Family | Acidibacter | NA |
| Seq96 | Bacteria | Proteobacteria | Gammaproteobacteria | R7C24 | NA | NA | NA |
| Seq170 | Bacteria | Proteobacteria | Gammaproteobacteria | Gammaproteobacteria Incertae Sedis | Unknown Family | Acidibacter | NA |
| Seq94 | Bacteria | Proteobacteria | Gammaproteobacteria | Xanthomonadales | Rhodanobacteraceae | Dokdonella | ginsengisoli |
| Seq114 | Bacteria | Proteobacteria | Gammaproteobacteria | Diplorickettsiales | Diplorickettsiaceae | NA | NA |
| Seq235 | Bacteria | Proteobacteria | Gammaproteobacteria | Legionellales | Legionellaceae | Legionella | NA |
$cyan \(cyan\)ps phyloseq-class experiment-level object otu_table() OTU Table: [ 30 taxa and 35 samples ] sample_data() Sample Data: [ 35 samples by 3 sample variables ] tax_table() Taxonomy Table: [ 30 taxa by 7 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 30 tips and 29 internal nodes ] refseq() DNAStringSet: [ 30 reference sequences ]
\(cyan\)amp ampvis2 object with 5 elements. Summary of OTU table: Samples OTUs Total#Reads Min#Reads Max#Reads Median#Reads 35 30 176073 336 15108 4911 Avg#Reads 5030.66
Assigned taxonomy: Kingdom Phylum Class Order Family Genus Species 30(100%) 30(100%) 30(100%) 30(100%) 30(100%) 29(96.67%) 2(6.67%)
Metadata variables: 4 SampleID, Day, Description, Group
\(cyan\)heat
\(cyan\)heat_rel
\(cyan\)tree
\(cyan\)taxa
| Kingdom | Phylum | Class | Order | Family | Genus | Species | |
|---|---|---|---|---|---|---|---|
| Seq98 | Bacteria | Proteobacteria | Gammaproteobacteria | Burkholderiales | Burkholderiaceae | Burkholderia | telluris |
| Seq64 | Bacteria | Proteobacteria | Alphaproteobacteria | Azospirillales | Inquilinaceae | Inquilinus | ginsengisoli |
| Seq4 | Bacteria | Proteobacteria | Alphaproteobacteria | Azospirillales | Inquilinaceae | Inquilinus | NA |
| Seq244 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhizobiales | Xanthobacteraceae | Pseudorhodoplanes | NA |
| Seq11 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhizobiales | Xanthobacteraceae | Bradyrhizobium | NA |
| Seq74 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhizobiales | Xanthobacteraceae | Bradyrhizobium | NA |
| Seq22 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhizobiales | Xanthobacteraceae | Starkeya | NA |
| Seq45 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhizobiales | Rhizobiaceae | Allorhizobium | NA |
| Seq9 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhizobiales | Rhizobiaceae | Allorhizobium | NA |
| Seq29 | Bacteria | Myxococcota | Polyangia | Nannocystales | Nannocystaceae | Nannocystis | NA |
| Seq239 | Bacteria | Myxococcota | Polyangia | Polyangiales | Polyangiaceae | Labilithrix | NA |
| Seq50 | Bacteria | Actinobacteriota | Actinobacteria | Corynebacteriales | Mycobacteriaceae | Mycobacterium | NA |
| Seq34 | Bacteria | Actinobacteriota | Actinobacteria | Streptosporangiales | Streptosporangiaceae | Herbidospora | NA |
| Seq56 | Bacteria | Firmicutes | Bacilli | Paenibacillales | Paenibacillaceae | Paenibacillus | NA |
| Seq25 | Bacteria | Firmicutes | Bacilli | Bacillales | Bacillaceae | Terribacillus | NA |
| Seq3 | Bacteria | Firmicutes | Bacilli | Bacillales | Bacillaceae | Bacillus | NA |
| Seq5 | Bacteria | Firmicutes | Bacilli | Bacillales | Bacillaceae | Bacillus | NA |
| Seq13 | Bacteria | Firmicutes | Bacilli | Bacillales | Planococcaceae | NA | NA |
| Seq7 | Bacteria | Firmicutes | Bacilli | Bacillales | Planococcaceae | Solibacillus | NA |
| Seq19 | Bacteria | Firmicutes | Bacilli | Bacillales | Planococcaceae | Solibacillus | NA |
| Seq17 | Bacteria | Planctomycetota | Planctomycetes | Isosphaerales | Isosphaeraceae | Singulisphaera | NA |
| Seq324 | Bacteria | Acidobacteriota | Acidobacteriae | Acidobacteriales | Acidobacteriaceae (Subgroup 1) | Edaphobacter | NA |
| Seq16 | Bacteria | Bacteroidota | Bacteroidia | Cytophagales | Microscillaceae | Ohtaekwangia | NA |
| Seq6 | Bacteria | Bacteroidota | Bacteroidia | Cytophagales | Microscillaceae | Ohtaekwangia | NA |
| Seq1 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Chitinophaga | NA |
| Seq109 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Chitinophaga | NA |
| Seq49 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Taibaiella | NA |
| Seq73 | Bacteria | Bacteroidota | Bacteroidia | Cytophagales | Spirosomaceae | Dyadobacter | NA |
| Seq15 | Bacteria | Proteobacteria | Gammaproteobacteria | Xanthomonadales | Rhodanobacteraceae | Luteibacter | NA |
| Seq26 | Bacteria | Proteobacteria | Gammaproteobacteria | Xanthomonadales | Xanthomonadaceae | Luteimonas | NA |
$darkgreen \(darkgreen\)ps phyloseq-class experiment-level object otu_table() OTU Table: [ 44 taxa and 35 samples ] sample_data() Sample Data: [ 35 samples by 3 sample variables ] tax_table() Taxonomy Table: [ 44 taxa by 7 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 44 tips and 43 internal nodes ] refseq() DNAStringSet: [ 44 reference sequences ]
\(darkgreen\)amp ampvis2 object with 5 elements. Summary of OTU table: Samples OTUs Total#Reads Min#Reads Max#Reads Median#Reads 35 44 12360 0 2468 65 Avg#Reads 353.14
Assigned taxonomy: Kingdom Phylum Class Order Family Genus Species 44(100%) 44(100%) 44(100%) 43(97.73%) 40(90.91%) 24(54.55%) 3(6.82%)
Metadata variables: 4 SampleID, Day, Description, Group
\(darkgreen\)heat
\(darkgreen\)heat_rel
\(darkgreen\)tree
\(darkgreen\)taxa
| Kingdom | Phylum | Class | Order | Family | Genus | Species | |
|---|---|---|---|---|---|---|---|
| Seq342 | Archaea | Crenarchaeota | Nitrososphaeria | Nitrososphaerales | Nitrososphaeraceae | NA | NA |
| Seq196 | Archaea | Crenarchaeota | Nitrososphaeria | Nitrososphaerales | Nitrososphaeraceae | Candidatus Nitrocosmicus | NA |
| Seq276 | Archaea | Crenarchaeota | Nitrososphaeria | Nitrososphaerales | Nitrososphaeraceae | NA | NA |
| Seq419 | Archaea | Crenarchaeota | Nitrososphaeria | Nitrososphaerales | Nitrososphaeraceae | NA | NA |
| Seq454 | Bacteria | Proteobacteria | Gammaproteobacteria | Diplorickettsiales | Diplorickettsiaceae | Aquicella | NA |
| Seq329 | Bacteria | Proteobacteria | Alphaproteobacteria | Sphingomonadales | Sphingomonadaceae | Sphingomonas | jaspsi |
| Seq332 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhodospirillales | Magnetospiraceae | NA | NA |
| Seq429 | Bacteria | Proteobacteria | Alphaproteobacteria | Acetobacterales | Acetobacteraceae | NA | NA |
| Seq195 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhizobiales | Xanthobacteraceae | Pseudorhodoplanes | NA |
| Seq144 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhizobiales | Xanthobacteraceae | Pseudolabrys | NA |
| Seq369 | Bacteria | Proteobacteria | Alphaproteobacteria | Micropepsales | Micropepsaceae | NA | NA |
| Seq362 | Bacteria | Proteobacteria | Alphaproteobacteria | Sphingomonadales | Sphingomonadaceae | Altererythrobacter | NA |
| Seq428 | Bacteria | Proteobacteria | Alphaproteobacteria | Caulobacterales | Caulobacteraceae | Phenylobacterium | NA |
| Seq347 | Bacteria | Nitrospirota | Nitrospiria | Nitrospirales | Nitrospiraceae | Nitrospira | japonica |
| Seq328 | Bacteria | Myxococcota | Polyangia | Haliangiales | Haliangiaceae | Haliangium | NA |
| Seq348 | Bacteria | Myxococcota | Polyangia | Polyangiales | BIrii41 | NA | NA |
| Seq308 | Bacteria | Myxococcota | Polyangia | Polyangiales | BIrii41 | NA | NA |
| Seq256 | Bacteria | Myxococcota | Polyangia | Polyangiales | Polyangiaceae | Minicystis | NA |
| Seq302 | Bacteria | Bdellovibrionota | Oligoflexia | 0319-6G20 | NA | NA | NA |
| Seq435 | Bacteria | Actinobacteriota | Actinobacteria | Propionibacteriales | Propionibacteriaceae | Microlunatus | NA |
| Seq291 | Bacteria | Actinobacteriota | Acidimicrobiia | IMCC26256 | NA | NA | NA |
| Seq392 | Bacteria | Actinobacteriota | Acidimicrobiia | Microtrichales | Ilumatobacteraceae | NA | NA |
| Seq321 | Bacteria | Bacteroidota | Bacteroidia | Sphingobacteriales | NS11-12 marine group | NA | NA |
| Seq186 | Bacteria | Firmicutes | Bacilli | Paenibacillales | Paenibacillaceae | Paenibacillus | NA |
| Seq204 | Bacteria | Firmicutes | Bacilli | Bacillales | Planococcaceae | Paenisporosarcina | NA |
| Seq146 | Bacteria | Spirochaetota | Spirochaetia | Spirochaetales | Spirochaetaceae | Spirochaeta 2 | NA |
| Seq359 | Bacteria | Patescibacteria | Saccharimonadia | Saccharimonadales | LWQ8 | NA | NA |
| Seq269 | Bacteria | Chloroflexi | Chloroflexia | Chloroflexales | Roseiflexaceae | NA | NA |
| Seq267 | Bacteria | Verrucomicrobiota | Verrucomicrobiae | Opitutales | Opitutaceae | Lacunisphaera | NA |
| Seq320 | Bacteria | Verrucomicrobiota | Verrucomicrobiae | Opitutales | Opitutaceae | Opitutus | NA |
| Seq203 | Bacteria | Verrucomicrobiota | Verrucomicrobiae | Chthoniobacterales | Terrimicrobiaceae | Terrimicrobium | NA |
| Seq497 | Bacteria | Planctomycetota | Planctomycetes | Pirellulales | Pirellulaceae | NA | NA |
| Seq384 | Bacteria | Planctomycetota | Planctomycetes | Planctomycetales | NA | NA | NA |
| Seq155 | Bacteria | Acidobacteriota | Blastocatellia | Blastocatellales | Blastocatellaceae | NA | NA |
| Seq350 | Bacteria | Bacteroidota | Bacteroidia | Sphingobacteriales | env.OPS 17 | NA | NA |
| Seq224 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Terrimonas | NA |
| Seq409 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Edaphobaculum | NA |
| Seq337 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Taibaiella | NA |
| Seq494 | Bacteria | Bacteroidota | Bacteroidia | Cytophagales | Amoebophilaceae | Candidatus Amoebophilus | NA |
| Seq277 | Bacteria | Proteobacteria | Gammaproteobacteria | NA | NA | NA | NA |
| Seq191 | Bacteria | Proteobacteria | Gammaproteobacteria | Gammaproteobacteria Incertae Sedis | Unknown Family | Acidibacter | NA |
| Seq335 | Bacteria | Proteobacteria | Gammaproteobacteria | Xanthomonadales | Xanthomonadaceae | Luteimonas | vadosa |
| Seq157 | Bacteria | Proteobacteria | Gammaproteobacteria | Xanthomonadales | Rhodanobacteraceae | Dokdonella | NA |
| Seq385 | Bacteria | Proteobacteria | Gammaproteobacteria | Salinisphaerales | Solimonadaceae | NA | NA |
$plum \(plum\)ps phyloseq-class experiment-level object otu_table() OTU Table: [ 30 taxa and 35 samples ] sample_data() Sample Data: [ 35 samples by 3 sample variables ] tax_table() Taxonomy Table: [ 30 taxa by 7 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 30 tips and 29 internal nodes ] refseq() DNAStringSet: [ 30 reference sequences ]
\(plum\)amp ampvis2 object with 5 elements. Summary of OTU table: Samples OTUs Total#Reads Min#Reads Max#Reads Median#Reads 35 30 17172 0 3034 103 Avg#Reads 490.63
Assigned taxonomy: Kingdom Phylum Class Order Family Genus Species 30(100%) 30(100%) 30(100%) 30(100%) 27(90%) 16(53.33%) 1(3.33%)
Metadata variables: 4 SampleID, Day, Description, Group
\(plum\)heat
\(plum\)heat_rel
\(plum\)tree
\(plum\)taxa
| Kingdom | Phylum | Class | Order | Family | Genus | Species | |
|---|---|---|---|---|---|---|---|
| Seq161 | Bacteria | Proteobacteria | Gammaproteobacteria | Diplorickettsiales | Diplorickettsiaceae | Aquicella | NA |
| Seq181 | Bacteria | Proteobacteria | Gammaproteobacteria | CCD24 | NA | NA | NA |
| Seq207 | Bacteria | Cyanobacteria | Vampirivibrionia | Obscuribacterales | Obscuribacteraceae | NA | NA |
| Seq219 | Bacteria | Proteobacteria | Alphaproteobacteria | Sphingomonadales | Sphingomonadaceae | Sphingomonas | naasensis |
| Seq216 | Bacteria | Proteobacteria | Alphaproteobacteria | Dongiales | Dongiaceae | Dongia | NA |
| Seq220 | Bacteria | Proteobacteria | Alphaproteobacteria | Reyranellales | Reyranellaceae | Reyranella | NA |
| Seq90 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhizobiales | Xanthobacteraceae | Pseudolabrys | NA |
| Seq271 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhizobiales | Xanthobacteraceae | NA | NA |
| Seq111 | Bacteria | Proteobacteria | Alphaproteobacteria | Sphingomonadales | Sphingomonadaceae | Sphingobium | NA |
| Seq33 | Bacteria | Proteobacteria | Alphaproteobacteria | Caulobacterales | Caulobacteraceae | Caulobacter | NA |
| Seq189 | Bacteria | Actinobacteriota | Thermoleophilia | Solirubrobacterales | Solirubrobacteraceae | Solirubrobacter | NA |
| Seq85 | Bacteria | Myxococcota | Polyangia | Polyangiales | BIrii41 | NA | NA |
| Seq208 | Bacteria | Actinobacteriota | Actinobacteria | Micrococcales | Microbacteriaceae | Galbitalea | NA |
| Seq218 | Bacteria | Bacteroidota | Bacteroidia | Sphingobacteriales | Sphingobacteriaceae | Mucilaginibacter | NA |
| Seq215 | Bacteria | Firmicutes | Bacilli | Bacillales | Bacillaceae | Bacillus | NA |
| Seq177 | Bacteria | Fibrobacterota | Fibrobacteria | Fibrobacterales | Fibrobacteraceae | NA | NA |
| Seq375 | Bacteria | Chloroflexi | Chloroflexia | Thermomicrobiales | JG30-KF-CM45 | NA | NA |
| Seq430 | Bacteria | Planctomycetota | Phycisphaerae | Phycisphaerales | Phycisphaeraceae | NA | NA |
| Seq368 | Bacteria | Planctomycetota | Planctomycetes | Pirellulales | Pirellulaceae | NA | NA |
| Seq152 | Bacteria | Planctomycetota | Planctomycetes | Gemmatales | Gemmataceae | Gemmata | NA |
| Seq304 | Bacteria | Planctomycetota | Planctomycetes | Planctomycetales | NA | NA | NA |
| Seq141 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Terrimonas | NA |
| Seq343 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | NA | NA |
| Seq412 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Edaphobaculum | NA |
| Seq457 | Bacteria | Proteobacteria | Gammaproteobacteria | Coxiellales | Coxiellaceae | Coxiella | NA |
| Seq130 | Bacteria | Proteobacteria | Gammaproteobacteria | R7C24 | NA | NA | NA |
| Seq346 | Bacteria | Proteobacteria | Gammaproteobacteria | Salinisphaerales | Solimonadaceae | Alkanibacter | NA |
| Seq474 | Bacteria | Proteobacteria | Gammaproteobacteria | Diplorickettsiales | Diplorickettsiaceae | NA | NA |
| Seq472 | Bacteria | Proteobacteria | Gammaproteobacteria | Diplorickettsiales | Diplorickettsiaceae | NA | NA |
| Seq227 | Bacteria | Proteobacteria | Gammaproteobacteria | Diplorickettsiales | Diplorickettsiaceae | NA | NA |
$red \(red\)ps phyloseq-class experiment-level object otu_table() OTU Table: [ 75 taxa and 35 samples ] sample_data() Sample Data: [ 35 samples by 3 sample variables ] tax_table() Taxonomy Table: [ 75 taxa by 7 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 75 tips and 74 internal nodes ] refseq() DNAStringSet: [ 75 reference sequences ]
\(red\)amp ampvis2 object with 5 elements. Summary of OTU table: Samples OTUs Total#Reads Min#Reads Max#Reads Median#Reads 35 75 66650 241 5339 1528 Avg#Reads 1904.29
Assigned taxonomy: Kingdom Phylum Class Order Family Genus Species 75(100%) 75(100%) 75(100%) 74(98.67%) 74(98.67%) 66(88%) 11(14.67%)
Metadata variables: 4 SampleID, Day, Description, Group
\(red\)heat
\(red\)heat_rel
\(red\)tree
\(red\)taxa
| Kingdom | Phylum | Class | Order | Family | Genus | Species | |
|---|---|---|---|---|---|---|---|
| Seq117 | Bacteria | Proteobacteria | Gammaproteobacteria | Burkholderiales | Burkholderiaceae | Burkholderia | NA |
| Seq60 | Bacteria | Proteobacteria | Gammaproteobacteria | Burkholderiales | Alcaligenaceae | NA | NA |
| Seq121 | Bacteria | Proteobacteria | Gammaproteobacteria | Burkholderiales | Comamonadaceae | Variovorax | NA |
| Seq383 | Bacteria | Proteobacteria | Alphaproteobacteria | Ferrovibrionales | Ferrovibrionaceae | Ferrovibrio | soli |
| Seq228 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhizobiales | Xanthobacteraceae | NA | NA |
| Seq178 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhizobiales | Beijerinckiaceae | Methylobacterium-Methylorubrum | NA |
| Seq171 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhizobiales | Beijerinckiaceae | Bosea | thiooxidans |
| Seq341 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhizobiales | Beijerinckiaceae | Bosea | NA |
| Seq89 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhizobiales | Xanthobacteraceae | Tardiphaga | robiniae |
| Seq31 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhizobiales | Xanthobacteraceae | Starkeya | NA |
| Seq160 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhizobiales | Hyphomicrobiaceae | Hyphomicrobium | NA |
| Seq251 | Bacteria | Proteobacteria | Alphaproteobacteria | Sphingomonadales | Sphingomonadaceae | Sphingopyxis | NA |
| Seq373 | Bacteria | Proteobacteria | Alphaproteobacteria | Sphingomonadales | Sphingomonadaceae | Altererythrobacter | NA |
| Seq108 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhizobiales | Rhizobiaceae | Mesorhizobium | NA |
| Seq71 | Bacteria | Actinobacteriota | Thermoleophilia | Solirubrobacterales | Solirubrobacteraceae | Conexibacter | NA |
| Seq338 | Bacteria | Myxococcota | Polyangia | Polyangiales | Polyangiaceae | Pajaroellobacter | NA |
| Seq230 | Bacteria | Myxococcota | Polyangia | Polyangiales | Polyangiaceae | Pajaroellobacter | NA |
| Seq88 | Bacteria | Myxococcota | Polyangia | Polyangiales | Polyangiaceae | Sorangium | NA |
| Seq303 | Bacteria | Actinobacteriota | Acidimicrobiia | Microtrichales | Ilumatobacteraceae | NA | NA |
| Seq93 | Bacteria | Actinobacteriota | Actinobacteria | Micrococcales | Promicromonosporaceae | Promicromonospora | NA |
| Seq366 | Bacteria | Actinobacteriota | Actinobacteria | Micrococcales | Microbacteriaceae | Leifsonia | NA |
| Seq134 | Bacteria | Actinobacteriota | Actinobacteria | Micrococcales | Microbacteriaceae | NA | aoyamense |
| Seq411 | Bacteria | Actinobacteriota | Actinobacteria | Micrococcales | Microbacteriaceae | Agromyces | NA |
| Seq115 | Bacteria | Actinobacteriota | Actinobacteria | Corynebacteriales | Mycobacteriaceae | Mycobacterium | NA |
| Seq128 | Bacteria | Actinobacteriota | Actinobacteria | Streptosporangiales | Streptosporangiaceae | Herbidospora | mongoliensis |
| Seq398 | Bacteria | Actinobacteriota | Actinobacteria | Streptosporangiales | Streptosporangiaceae | Nonomuraea | NA |
| Seq78 | Bacteria | Actinobacteriota | Actinobacteria | Streptosporangiales | Thermomonosporaceae | Actinocorallia | NA |
| Seq193 | Bacteria | Actinobacteriota | Actinobacteria | Streptomycetales | Streptomycetaceae | Streptomyces | NA |
| Seq148 | Bacteria | Bacteroidota | Bacteroidia | Sphingobacteriales | Sphingobacteriaceae | Pedobacter | NA |
| Seq307 | Bacteria | Bacteroidota | Bacteroidia | Sphingobacteriales | NS11-12 marine group | NA | NA |
| Seq431 | Bacteria | Bacteroidota | Bacteroidia | Sphingobacteriales | NS11-12 marine group | NA | NA |
| Seq205 | Bacteria | Firmicutes | Bacilli | Paenibacillales | Paenibacillaceae | Paenibacillus | lautus |
| Seq86 | Bacteria | Firmicutes | Bacilli | Paenibacillales | Paenibacillaceae | Paenibacillus | NA |
| Seq51 | Bacteria | Firmicutes | Bacilli | Paenibacillales | Paenibacillaceae | Paenibacillus | NA |
| Seq18 | Bacteria | Firmicutes | Bacilli | Paenibacillales | Paenibacillaceae | Paenibacillus | NA |
| Seq192 | Bacteria | Firmicutes | Bacilli | Paenibacillales | Paenibacillaceae | Paenibacillus | NA |
| Seq266 | Bacteria | Firmicutes | Bacilli | Paenibacillales | Paenibacillaceae | Paenibacillus | NA |
| Seq278 | Bacteria | Firmicutes | Bacilli | Bacillales | Planococcaceae | Domibacillus | NA |
| Seq99 | Bacteria | Firmicutes | Bacilli | Bacillales | Bacillaceae | Bacillus | NA |
| Seq53 | Bacteria | Firmicutes | Bacilli | Bacillales | Planococcaceae | Lysinibacillus | NA |
| Seq118 | Bacteria | Firmicutes | Bacilli | Bacillales | Planococcaceae | Lysinibacillus | NA |
| Seq76 | Bacteria | Verrucomicrobiota | Verrucomicrobiae | Chthoniobacterales | Terrimicrobiaceae | Terrimicrobium | NA |
| Seq232 | Bacteria | Verrucomicrobiota | Verrucomicrobiae | Chthoniobacterales | Terrimicrobiaceae | Terrimicrobium | NA |
| Seq168 | Bacteria | Verrucomicrobiota | Verrucomicrobiae | Chthoniobacterales | Terrimicrobiaceae | Terrimicrobium | NA |
| Seq387 | Bacteria | Verrucomicrobiota | Verrucomicrobiae | Chthoniobacterales | Chthoniobacteraceae | Chthoniobacter | NA |
| Seq285 | Bacteria | Verrucomicrobiota | Verrucomicrobiae | Chthoniobacterales | Chthoniobacteraceae | Chthoniobacter | NA |
| Seq217 | Bacteria | Verrucomicrobiota | Verrucomicrobiae | Verrucomicrobiales | Verrucomicrobiaceae | Verrucomicrobium | spinosum |
| Seq190 | Bacteria | Planctomycetota | Planctomycetes | Gemmatales | Gemmataceae | Gemmata | NA |
| Seq249 | Bacteria | Planctomycetota | Planctomycetes | Isosphaerales | Isosphaeraceae | Singulisphaera | NA |
| Seq254 | Bacteria | Acidobacteriota | Acidobacteriae | Acidobacteriales | Acidobacteriaceae (Subgroup 1) | Edaphobacter | NA |
| Seq67 | Bacteria | Bacteroidota | Bacteroidia | Cytophagales | Microscillaceae | Ohtaekwangia | NA |
| Seq176 | Bacteria | Bacteroidota | Bacteroidia | Cytophagales | Microscillaceae | Ohtaekwangia | NA |
| Seq268 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Pseudoflavitalea | NA |
| Seq32 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Pseudoflavitalea | NA |
| Seq84 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Pseudoflavitalea | NA |
| Seq459 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Niastella | NA |
| Seq40 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Niastella | hibisci |
| Seq173 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Niastella | NA |
| Seq246 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Niastella | NA |
| Seq112 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Flavitalea | NA |
| Seq126 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Chitinophaga | ginsengihumi |
| Seq102 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Chitinophaga | arvensicola |
| Seq309 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Chitinophaga | soli |
| Seq70 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Chitinophaga | NA |
| Seq262 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Chitinophaga | NA |
| Seq301 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Taibaiella | NA |
| Seq572 | Bacteria | Proteobacteria | Gammaproteobacteria | Pseudomonadales | Moraxellaceae | NA | NA |
| Seq436 | Bacteria | Proteobacteria | Gammaproteobacteria | Steroidobacterales | Steroidobacteraceae | Steroidobacter | NA |
| Seq261 | Bacteria | Proteobacteria | Gammaproteobacteria | Gammaproteobacteria Incertae Sedis | Unknown Family | Acidibacter | NA |
| Seq298 | Bacteria | Proteobacteria | Gammaproteobacteria | Gammaproteobacteria Incertae Sedis | Unknown Family | Acidibacter | NA |
| Seq354 | Bacteria | Proteobacteria | Gammaproteobacteria | Gammaproteobacteria Incertae Sedis | Unknown Family | Acidibacter | NA |
| Seq226 | Bacteria | Proteobacteria | Gammaproteobacteria | Xanthomonadales | Rhodanobacteraceae | Tahibacter | NA |
| Seq147 | Bacteria | Proteobacteria | Gammaproteobacteria | Xanthomonadales | Xanthomonadaceae | Xanthomonas | NA |
| Seq370 | Bacteria | Proteobacteria | Gammaproteobacteria | Diplorickettsiales | Diplorickettsiaceae | NA | NA |
| Seq486 | Bacteria | Proteobacteria | Gammaproteobacteria | NA | NA | NA | NA |
$salmon \(salmon\)ps phyloseq-class experiment-level object otu_table() OTU Table: [ 73 taxa and 35 samples ] sample_data() Sample Data: [ 35 samples by 3 sample variables ] tax_table() Taxonomy Table: [ 73 taxa by 7 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 73 tips and 72 internal nodes ] refseq() DNAStringSet: [ 73 reference sequences ]
\(salmon\)amp ampvis2 object with 5 elements. Summary of OTU table: Samples OTUs Total#Reads Min#Reads Max#Reads Median#Reads 35 73 115332 0 20534 779 Avg#Reads 3295.2
Assigned taxonomy: Kingdom Phylum Class Order Family Genus Species 73(100%) 73(100%) 73(100%) 72(98.63%) 71(97.26%) 66(90.41%) 12(16.44%)
Metadata variables: 4 SampleID, Day, Description, Group
\(salmon\)heat
\(salmon\)heat_rel
\(salmon\)tree
\(salmon\)taxa
| Kingdom | Phylum | Class | Order | Family | Genus | Species | |
|---|---|---|---|---|---|---|---|
| Seq79 | Bacteria | Proteobacteria | Gammaproteobacteria | Burkholderiales | Oxalobacteraceae | Massilia | armeniaca |
| Seq122 | Bacteria | Proteobacteria | Gammaproteobacteria | Burkholderiales | Oxalobacteraceae | NA | NA |
| Seq179 | Bacteria | Proteobacteria | Gammaproteobacteria | Burkholderiales | Oxalobacteraceae | Pseudoduganella | NA |
| Seq145 | Bacteria | Proteobacteria | Gammaproteobacteria | Burkholderiales | Oxalobacteraceae | Pseudoduganella | eburnea |
| Seq154 | Bacteria | Proteobacteria | Gammaproteobacteria | Burkholderiales | Oxalobacteraceae | Massilia | NA |
| Seq222 | Bacteria | Proteobacteria | Gammaproteobacteria | Burkholderiales | Oxalobacteraceae | Massilia | NA |
| Seq159 | Bacteria | Proteobacteria | Gammaproteobacteria | Burkholderiales | Burkholderiaceae | Cupriavidus | NA |
| Seq8 | Bacteria | Proteobacteria | Gammaproteobacteria | Burkholderiales | Burkholderiaceae | Cupriavidus | NA |
| Seq14 | Bacteria | Proteobacteria | Gammaproteobacteria | Burkholderiales | Alcaligenaceae | Achromobacter | NA |
| Seq54 | Bacteria | Proteobacteria | Gammaproteobacteria | Burkholderiales | Comamonadaceae | NA | paradoxus |
| Seq95 | Bacteria | Proteobacteria | Gammaproteobacteria | Burkholderiales | Comamonadaceae | Xylophilus | NA |
| Seq340 | Bacteria | Proteobacteria | Gammaproteobacteria | Burkholderiales | Comamonadaceae | Rhizobacter | NA |
| Seq264 | Bacteria | Proteobacteria | Alphaproteobacteria | Sphingomonadales | Sphingomonadaceae | Sphingomonas | mucosissima |
| Seq333 | Bacteria | Proteobacteria | Alphaproteobacteria | Reyranellales | Reyranellaceae | Reyranella | NA |
| Seq158 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhizobiales | Beijerinckiaceae | Microvirga | NA |
| Seq248 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhizobiales | Rhizobiaceae | Allorhizobium | NA |
| Seq245 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhizobiales | Devosiaceae | Devosia | neptuniae |
| Seq185 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhizobiales | Rhizobiaceae | Shinella | NA |
| Seq156 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhizobiales | Rhizobiaceae | NA | NA |
| Seq42 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhizobiales | Rhizobiaceae | NA | NA |
| Seq286 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhizobiales | Rhizobiaceae | Ensifer | NA |
| Seq21 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhizobiales | Rhizobiaceae | Allorhizobium | NA |
| Seq107 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhizobiales | Rhizobiaceae | Neorhizobium | NA |
| Seq231 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhizobiales | Rhizobiaceae | Allorhizobium | azooxidifex |
| Seq140 | Bacteria | Proteobacteria | Alphaproteobacteria | Rhizobiales | Rhizobiaceae | NA | NA |
| Seq378 | Bacteria | Proteobacteria | Alphaproteobacteria | Caulobacterales | Caulobacteraceae | Phenylobacterium | mobile |
| Seq374 | Bacteria | Myxococcota | Polyangia | Haliangiales | Haliangiaceae | Haliangium | NA |
| Seq257 | Bacteria | Myxococcota | Polyangia | Polyangiales | Polyangiaceae | Pajaroellobacter | NA |
| Seq223 | Bacteria | Actinobacteriota | Actinobacteria | Micrococcales | Promicromonosporaceae | Cellulosimicrobium | NA |
| Seq83 | Bacteria | Actinobacteriota | Actinobacteria | Micrococcales | Microbacteriaceae | Microbacterium | NA |
| Seq68 | Bacteria | Actinobacteriota | Actinobacteria | Corynebacteriales | Mycobacteriaceae | Mycobacterium | NA |
| Seq293 | Bacteria | Actinobacteriota | Actinobacteria | Glycomycetales | Glycomycetaceae | Glycomyces | NA |
| Seq72 | Bacteria | Actinobacteriota | Actinobacteria | Streptomycetales | Streptomycetaceae | Streptomyces | NA |
| Seq101 | Bacteria | Actinobacteriota | Actinobacteria | Streptomycetales | Streptomycetaceae | Streptomyces | NA |
| Seq265 | Bacteria | Actinobacteriota | Actinobacteria | Streptomycetales | Streptomycetaceae | Streptomyces | NA |
| Seq260 | Bacteria | Bacteroidota | Bacteroidia | Sphingobacteriales | NA | NA | NA |
| Seq180 | Bacteria | Bacteroidota | Bacteroidia | Sphingobacteriales | Sphingobacteriaceae | Pedobacter | panaciterrae |
| Seq142 | Bacteria | Verrucomicrobiota | Verrucomicrobiae | Verrucomicrobiales | Verrucomicrobiaceae | Roseimicrobium | NA |
| Seq163 | Bacteria | Bacteroidota | Bacteroidia | Cytophagales | Microscillaceae | Ohtaekwangia | NA |
| Seq253 | Bacteria | Bacteroidota | Bacteroidia | Cytophagales | Microscillaceae | Ohtaekwangia | NA |
| Seq116 | Bacteria | Bacteroidota | Bacteroidia | Cytophagales | Microscillaceae | Ohtaekwangia | NA |
| Seq184 | Bacteria | Bacteroidota | Bacteroidia | Cytophagales | Microscillaceae | Ohtaekwangia | NA |
| Seq52 | Bacteria | Bacteroidota | Bacteroidia | Flavobacteriales | Flavobacteriaceae | Flavobacterium | NA |
| Seq187 | Bacteria | Bacteroidota | Bacteroidia | Flavobacteriales | Weeksellaceae | Chryseobacterium | ginsenosidimutans |
| Seq172 | Bacteria | Bacteroidota | Bacteroidia | Flavobacteriales | Weeksellaceae | Chryseobacterium | NA |
| Seq283 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Flavitalea | NA |
| Seq113 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Niastella | NA |
| Seq110 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Niastella | NA |
| Seq132 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Pseudoflavitalea | NA |
| Seq37 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Chitinophaga | NA |
| Seq162 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Chitinophaga | NA |
| Seq182 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Chitinophaga | humicola |
| Seq92 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Chitinophaga | NA |
| Seq80 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Chitinophaga | pinensis |
| Seq2 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Chitinophaga | NA |
| Seq221 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Chitinophaga | NA |
| Seq20 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Chitinophaga | NA |
| Seq282 | Bacteria | Bacteroidota | Bacteroidia | Chitinophagales | Chitinophagaceae | Chitinophaga | NA |
| Seq38 | Bacteria | Bacteroidota | Bacteroidia | Cytophagales | Spirosomaceae | Dyadobacter | NA |
| Seq237 | Bacteria | Proteobacteria | Gammaproteobacteria | Steroidobacterales | Steroidobacteraceae | Steroidobacter | agariperforans |
| Seq234 | Bacteria | Proteobacteria | Gammaproteobacteria | Xanthomonadales | Xanthomonadaceae | Pseudoxanthomonas | NA |
| Seq59 | Bacteria | Proteobacteria | Gammaproteobacteria | Xanthomonadales | Xanthomonadaceae | Stenotrophomonas | NA |
| Seq105 | Bacteria | Proteobacteria | Gammaproteobacteria | Xanthomonadales | Xanthomonadaceae | Stenotrophomonas | NA |
| Seq344 | Bacteria | Proteobacteria | Gammaproteobacteria | Xanthomonadales | Xanthomonadaceae | Lysobacter | NA |
| Seq41 | Bacteria | Proteobacteria | Gammaproteobacteria | Xanthomonadales | Xanthomonadaceae | Lysobacter | NA |
| Seq243 | Bacteria | Proteobacteria | Gammaproteobacteria | NA | NA | NA | NA |
| Seq406 | Bacteria | Proteobacteria | Gammaproteobacteria | Legionellales | Legionellaceae | Legionella | NA |
| Seq139 | Bacteria | Proteobacteria | Gammaproteobacteria | Pseudomonadales | Pseudomonadaceae | Pseudomonas | NA |
| Seq24 | Bacteria | Proteobacteria | Gammaproteobacteria | Pseudomonadales | Pseudomonadaceae | Pseudomonas | NA |
| Seq27 | Bacteria | Proteobacteria | Gammaproteobacteria | Pseudomonadales | Pseudomonadaceae | Pseudomonas | NA |
| Seq87 | Bacteria | Proteobacteria | Gammaproteobacteria | Pseudomonadales | Pseudomonadaceae | Pseudomonas | NA |
| Seq55 | Bacteria | Proteobacteria | Gammaproteobacteria | Pseudomonadales | Pseudomonadaceae | Pseudomonas | NA |
| Seq201 | Bacteria | Proteobacteria | Gammaproteobacteria | Enterobacterales | Enterobacteriaceae | Klebsiella | NA |
mpd по стадиям
physeq_merged <- merge_samples(ps.f, "Group", fun=sum)
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
ps.f@sam_data
## Day Description Group
## straw-16s-D01-1 D01 D01_1 early
## straw-16s-D01-2 D01 D01_2 early
## straw-16s-D01-3 D01 D01_3 early
## straw-16s-D03-1 D03 D03_1 early
## straw-16s-D03-2 D03 D03_2 early
## straw-16s-D03-3 D03 D03_3 early
## straw-16s-D05-1 D05 D05_1 early
## straw-16s-D05-2 D05 D05_2 early
## straw-16s-D05-3 D05 D05_3 early
## straw-16s-D05-4 D05 D05_4 early
## straw-16s-D05-5 D05 D05_5 early
## straw-16s-D07-2 D07 D07_2 middle
## straw-16s-D07-3 D07 D07_3 middle
## straw-16s-D08-1 D08 D08_1 middle
## straw-16s-D08-2 D08 D08_2 middle
## straw-16s-D08-3 D08 D08_3 middle
## straw-16s-D10-1 D10 D10_1 middle
## straw-16s-D10-2 D10 D10_2 middle
## straw-16s-D10-3 D10 D10_3 middle
## straw-16s-D10-4 D10 D10_4 middle
## straw-16s-D10-5 D10 D10_5 middle
## straw-16s-D12-1 D12 D12_1 late
## straw-16s-D12-2 D12 D12_2 late
## straw-16s-D12-3 D12 D12_3 late
## straw-16s-D13-1 D13 D13_1 late
## straw-16s-D13-2 D13 D13_2 late
## straw-16s-D13-3 D13 D13_3 late
## straw-16s-D14-1 D14 D14_1 late
## straw-16s-D14-2 D14 D14_2 late
## straw-16s-D14-3 D14 D14_3 late
## straw-16s-D14-4 D14 D14_4 late
## straw-16s-D14-5 D14 D14_5 late
## straw-16s-D15-1 D15 D15_1 late
## straw-16s-D15-2 D15 D15_2 late
## straw-16s-D15-3 D15 D15_3 late
# picante::mpd(l_vst$blue$ps@otu_table@.Data, cophenetic(l_vst$blue$ps@phy_tree)) %>%
# mean(na.rm = TRUE)
#
# picante::mpd(l_vst$salmon$ps@otu_table@.Data, cophenetic(l_vst$salmon$ps@phy_tree)) %>%
# mean(na.rm = TRUE)
picante::ses.mpd(physeq_merged@otu_table@.Data, cophenetic(physeq_merged@phy_tree))
## ntaxa mpd.obs mpd.rand.mean mpd.rand.sd mpd.obs.rank mpd.obs.z
## early 364 1.742096 1.862383 0.02930118 1 -4.1051776
## middle 490 1.839665 1.863938 0.02298162 147 -1.0561720
## late 822 1.877129 1.863599 0.01354403 842 0.9989395
## mpd.obs.p runs
## early 0.001 999
## middle 0.147 999
## late 0.842 999
list.files("meta/")
## [1] "cell_realtime_stat.xlsx" "cell_resp_ch_stat.xlsx"
## [3] "period_legend.xlsx"
period_legend - соответствие номеров мешочков в 16S с днями, думаю лучше заменить везде обозначения 1, 3 и тд на дни cell_resp - данные по дыханию, по ним какую-нибудь простую статистику. Да, повторностей для эксперимента и контроля разное количество cell_realtime - это циклы выхода целлюлаз по реалтайму, я их нормировала по 16S, есть в этой же таблице. Что с ними делать особо не знаю, вроде вообще решили выкинуть
realtime.data <- readxl::read_excel("meta/cell_realtime_stat.xlsx")
period.data <- readxl::read_excel("meta/period_legend.xlsx")
resp.data <- readxl::read_excel("meta/cell_resp_ch_stat.xlsx")
realtime.data
## # A tibble: 36 × 29
## id day contr…¹ CELL_…² CELL_…³ CELL_…⁴ CELL_…⁵ CELL_…⁶ CELL_…⁷ CELL_…⁸
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 G01-1-1 3 15.6 33.2 31.1 33.6 36.6 33.8 38.0 33.0
## 2 G01-1-2 3 15.5 33.3 31.3 32.9 36.2 33.8 41 32.7
## 3 G01-1-3 3 15.5 33.6 31.2 32.9 35.8 33.3 38.2 32.9
## 4 G01-2-1 3 17.1 36.0 31.1 33.1 36.9 36.4 38.8 33.8
## 5 G01-2-2 3 17.3 35.6 31.0 32.8 37.8 35.7 40.1 33.6
## 6 G01-2-3 3 17.1 36 31.1 32.8 38.1 35.2 37.6 33.6
## 7 G01-3-1 3 16.4 35.4 31.6 34.9 39.1 35.3 40.9 33.2
## 8 G01-3-2 3 16.5 35.6 31.3 34.4 37.9 34.6 38.9 33.1
## 9 G01-3-3 3 16.4 35.8 31.5 34.4 41.6 35.0 36.5 33.4
## 10 G05-1-1 28 18.2 26.3 29.1 29.4 28.2 28.4 37.4 34.5
## # … with 26 more rows, 19 more variables: CELL_193122 <dbl>, CELL_73229 <dbl>,
## # CELL_47814 <dbl>, CELL_163125 <dbl>, CELL_73266 <dbl>, CELL_88582 <dbl>,
## # CELL_63583 <dbl>, CELL_14199 <dbl>, CELL_95616 <dbl>, CELL_63504 <dbl>,
## # CELL_08643 <chr>, CELL_199599 <dbl>, CELL_01426 <dbl>, CELL_71601 <dbl>,
## # CELL_45099 <dbl>, CELL_191900 <dbl>, CELL_99463 <dbl>, CELL_74579 <dbl>,
## # CELL_183403 <dbl>, and abbreviated variable names ¹contr_16S, ²CELL_172283,
## # ³CELL_203163, ⁴CELL_83325, ⁵CELL_188413, ⁶CELL_109631, ⁷CELL_188119, …
impute.mean <- function(x) replace(x, is.na(x), mean(x, na.rm = TRUE))
realtime_data <- realtime.data %>%
mutate(CELL_08643 = as.numeric(CELL_08643)) %>%
group_by(day) %>%
mutate(nice_cell = impute.mean(CELL_08643)) %>%
mutate(day = as.factor(day))
## Warning in mask$eval_all_mutate(quo): NAs introduced by coercion
realtime_zjena <- readxl::read_excel("cellulases_gene_expression(1).xlsx")
geom_mean <-function(x){exp(mean(log(x)))}
realtime_zjena_geom <- realtime_zjena %>%
mutate(repeats = paste0(realtime_zjena$week, "-", rep(1:3, 3, each=3))) %>%
relocate(repeats, 1) %>%
group_by(repeats) %>%
summarise_if(is.numeric, geom_mean) %>%
mutate(week = as.factor(week)) %>%
arrange(week)
UPGMA целлюлаз(дистанция - брей).
Не очень понимаю, что это значит. В зависимости от дистанции получаются
совсем разные кластеризации.
realtime_matrix <- realtime_zjena_geom %>%
column_to_rownames("repeats") %>%
select_if(is.numeric) %>%
as.matrix()
hcl <- hclust(vegan::vegdist(t(realtime_matrix), method="bray"), "average")
plot(hcl)
А вот - евклидова дистанция.
Ничего общего с предыдущей картинкой
hcl <- hclust(vegan::vegdist(t(realtime_matrix), method="euclidian"), "average")
plot(hcl)
Корплот целлюлаз.
Только корреляции
m = cor(realtime_matrix)
corrplot::corrplot(m)
Только достоверные корреляции
В общем - кластеры есть - но корреляцие недостоверные(за некоторыми
исключениями)
cor_test_mat <- psych::corr.test(realtime_matrix)$p
corrplot::corrplot(m, p.mat = cor_test_mat, method = 'circle', type = 'lower', insig='blank',
order = 'AOE', diag = FALSE)$corrPos -> p1
# text(p1$x, p1$y, round(p1$corr, 2))
То же, но матрица уже прологарифмирована.
cor_test_mat <- psych::corr.test(log(realtime_matrix))
corrplot::corrplot(cor_test_mat$r, p.mat = cor_test_mat$p, method = 'circle', type = 'lower', insig='blank',
order = 'AOE', diag = FALSE)$corrPos -> p1
# text(p1$x, p1$y, round(p1$corr, 2))
add phylogenic tree (mafft - iqtree(ML))
realtime_tree <- ape::read.tree("al_chz.fasta.contree")
plot(realtime_tree)
library(tidyverse)
realtime_data %>%
select(c("day", "id", "contr_16S", "CELL_172283")) %>%
mutate(bio_repl = gsub("-[1-3]$", "", id)) %>%
group_by(bio_repl, day) %>%
summarise(contr_tmean = mean(contr_16S),
data_tmean = mean(CELL_172283)) %>%
mutate(dCt = data_tmean - contr_tmean)
## `summarise()` has grouped output by 'bio_repl'. You can override using the
## `.groups` argument.
## # A tibble: 12 × 5
## # Groups: bio_repl [12]
## bio_repl day contr_tmean data_tmean dCt
## <chr> <fct> <dbl> <dbl> <dbl>
## 1 G01-1 3 15.5 33.4 17.8
## 2 G01-2 3 17.2 35.9 18.7
## 3 G01-3 3 16.4 35.6 19.2
## 4 G05-1 28 18.4 26.2 7.81
## 5 G05-2 28 16.5 26.2 9.73
## 6 G05-3 28 18.5 32.0 13.6
## 7 G10-1 91 17.4 26.2 8.72
## 8 G10-2 91 16.3 24.1 7.79
## 9 G10-3 91 17.5 26.3 8.81
## 10 G14-1 161 18.0 28.1 10.1
## 11 G14-2 161 16.5 27.0 10.6
## 12 G14-3 161 16.9 25.4 8.41
resp_data <- resp.data %>%
group_by(day) %>%
mutate(control = impute.mean(control)) %>%
mutate(straw = impute.mean(straw))
resp_data
## # A tibble: 164 × 3
## # Groups: day [27]
## day control straw
## <dbl> <dbl> <dbl>
## 1 0 250 250
## 2 3 285. 723.
## 3 3 263. 832.
## 4 3 131. 657.
## 5 3 307. 525.
## 6 3 241. 744.
## 7 3 245. 701.
## 8 3 245. 657.
## 9 7 274. 690.
## 10 7 296. 690.
## # … with 154 more rows
Дыхание - median(straw)/median(control)
Mожно использовать сторонний пакет(KneeArrower) что понять в каком месте knee_plot происходит перелом. Для элиминации повторностей возьмем медиану. Я хз как этот пакет работает(ищет производную, но как сглаживает - хз, там матан), но он говорит что перелом происходит скорее на 60-80 днях.
(я исправил ошибки - стало ближе к твоим данным)
https://github.com/agentlans/KneeArrower - вот здесь можно почитать про матан
# resp_data %>%
# filter(!day == 0) %>%
# group_by(day) %>%
# summarise(median_control = median(control),
# median_straw = median(straw)) %>%
# mutate(rel = median_straw/median_control) %>%
# ggplot() +
# geom_point(aes(x = day, y = rel))
resp_median <- resp_data %>%
filter(!day == 0) %>%
group_by(day) %>%
summarise(median_control = median(control),
median_straw = median(straw)) %>%
mutate(rel = median_straw/median_control)
thresholds <- c(0.25, 0.5, 0.75, 1)
# Find cutoff points at each threshold
cutoff.points <- lapply(thresholds, function(i) {
findCutoff(resp_median$day, resp_median$rel, method="first", i)
})
x.coord <- sapply(cutoff.points, function(p) p$x)
y.coord <- sapply(cutoff.points, function(p) p$y)
# Plot the cutoff points on the scatterplot
plot(resp_median$day, resp_median$rel, pch=20, col="gray")
points(x.coord, y.coord, col="red", pch=20)
text(x.coord, y.coord, labels=thresholds, pos=4, col="red")
period_data <- period.data %>%
mutate(bag_id = as.factor(bag_id))
period_data
## # A tibble: 15 × 2
## bag_id day
## <fct> <dbl>
## 1 1 3
## 2 2 7
## 3 3 14
## 4 4 21
## 5 5 28
## 6 6 35
## 7 7 49
## 8 8 63
## 9 9 77
## 10 10 91
## 11 11 105
## 12 12 119
## 13 13 140
## 14 14 161
## 15 15 182
Ну я вот не очень понимаю что делать дальше - привязать кластеры к этой картинке?
resp_median_bags <- resp_median %>%
left_join(period.data, by="day") %>%
mutate(bag_id = as.factor(bag_id))
ps.f.r <- rarefy_even_depth(ps.f)
## You set `rngseed` to FALSE. Make sure you've set & recorded
## the random seed of your session for reproducibility.
## See `?set.seed`
## ...
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
## 94OTUs were removed because they are no longer
## present in any sample after random subsampling
## ...
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
estimate_richness(ps.f.r) %>%
rownames_to_column("ID") %>%
mutate(bag_id = as.factor(
as.numeric(
gsub("\\..+","",
gsub("straw\\.16s\\.D","", ID)
)
)
)
) %>%
group_by(bag_id) %>%
summarise(Observed = mean(Observed),
Shannon = mean(Shannon),
InvSimpson = mean(InvSimpson)) %>%
left_join(period_data, by="bag_id") %>%
left_join(resp_median_bags, by="bag_id") %>%
mutate(
Observed_scaled = scale(Observed),
Shannon_scaled = scale(Shannon),
InvSimpson_scaled = scale(InvSimpson),
Respiration_scaled = scale(rel)
) %>%
select(c(bag_id, Observed_scaled, Shannon_scaled, InvSimpson_scaled, Respiration_scaled)) %>%
pivot_longer(c("Observed_scaled", "Shannon_scaled", "InvSimpson_scaled", "Respiration_scaled")) %>%
ggplot(aes(y = value, x = bag_id, group = name)) +
geom_line(aes(color = name),
alpha = 0.8) +
theme_bw()
Корреляция -ртрицательная слабенькая, не особо достоверная и только Пирсона(которая работает для нормального распределения(у нас хроносерия, надо спирмана по хорошему))
alpha_resp <- phyloseq::estimate_richness(ps.f.r) %>%
rownames_to_column("ID") %>%
mutate(bag_id = as.factor(
as.numeric(
gsub("\\..+","",
gsub("straw\\.16s\\.D","", ID)
)
)
)
) %>%
group_by(bag_id) %>%
summarise(Observed = mean(Observed),
Shannon = mean(Shannon),
InvSimpson = mean(InvSimpson)) %>%
left_join(period_data, by="bag_id") %>%
left_join(resp_median_bags, by="bag_id") %>%
mutate(
Observed_scaled = scale(Observed),
Shannon_scaled = scale(Shannon),
InvSimpson_scaled = scale(InvSimpson),
Respiration_scaled = scale(rel)
) %>%
select(c(bag_id, Observed_scaled, Shannon_scaled, InvSimpson_scaled, Respiration_scaled))
cor.test(alpha_resp$Observed_scaled, alpha_resp$Respiration_scaled, method = "spearman")
##
## Spearman's rank correlation rho
##
## data: alpha_resp$Observed_scaled and alpha_resp$Respiration_scaled
## S = 230, p-value = 0.2629
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.3939394
cor.test(alpha_resp$Shannon_scaled, alpha_resp$Respiration_scaled, method = "spearman")
##
## Spearman's rank correlation rho
##
## data: alpha_resp$Shannon_scaled and alpha_resp$Respiration_scaled
## S = 232, p-value = 0.2474
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.4060606
cor.test(alpha_resp$InvSimpson_scaled, alpha_resp$Respiration_scaled, method = "spearman")
##
## Spearman's rank correlation rho
##
## data: alpha_resp$InvSimpson_scaled and alpha_resp$Respiration_scaled
## S = 226, p-value = 0.2956
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.369697
cor.test(alpha_resp$Observed_scaled, alpha_resp$Respiration_scaled, method = "pearson")
##
## Pearson's product-moment correlation
##
## data: alpha_resp$Observed_scaled and alpha_resp$Respiration_scaled
## t = -2.7571, df = 8, p-value = 0.02479
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.9223009 -0.1220128
## sample estimates:
## cor
## -0.6980158
cor.test(alpha_resp$Shannon_scaled, alpha_resp$Respiration_scaled, method = "pearson")
##
## Pearson's product-moment correlation
##
## data: alpha_resp$Shannon_scaled and alpha_resp$Respiration_scaled
## t = -2.8453, df = 8, p-value = 0.02164
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.9255437 -0.1437787
## sample estimates:
## cor
## -0.7092031
cor.test(alpha_resp$InvSimpson_scaled, alpha_resp$Respiration_scaled, method = "pearson")
##
## Pearson's product-moment correlation
##
## data: alpha_resp$InvSimpson_scaled and alpha_resp$Respiration_scaled
## t = -2.2565, df = 8, p-value = 0.05401
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.899910986 0.009841998
## sample estimates:
## cor
## -0.6236489
В табличке есть колонка cluster - “exl” в ней обозначает что это та часть датасета которая не пошла в WGCNA - мажоры в единичных семплах.
ps.m <- phyloseq::psmelt(ps.f)
ps.m <- ps.m %>%
mutate_if(is.character, as.factor)
ps.data.out <- ps.m %>%
select(-Group) %>%
pivot_wider(names_from = c(Day, Description, Sample), values_from = Abundance, values_fill = 0)
#create empty dataframe with columnnames
external_empty_dataframe <- data.frame(OTU=factor(), cluster=factor(), stringsAsFactors = FALSE)
for (i in names(l_vst)) {
a <- taxa_names(l_vst[[i]][["ps"]])
b <- rep(i, length(a))
d <- data.frame(OTU = as.factor(a),
cluster = as.factor(b))
external_empty_dataframe <- rbind(external_empty_dataframe, d)
}
clusters.otu.df <- external_empty_dataframe
# add exl taxa -- taxa "exl"
d <- data.frame(OTU = as.factor(ps.exl.taxa),
cluster = as.factor(rep("exl", length(ps.exl.taxa))))
clusters.otu.df <- rbind(clusters.otu.df, d)
ps.data.out.exl <- left_join(clusters.otu.df, ps.data.out, by="OTU")
# write.table(ps.data.out.exl, file = "ps.data.out.tsv", sep = "\t")
ps.data.out.exl
Bdellovibrionota - хищники, индикатор развитого сообщества
ps.m %>%
filter(Phylum == "Bdellovibrionota") %>%
group_by(Description, Day) %>%
summarise(Bs = sum(Abundance)) %>%
ggplot() +
geom_boxplot(aes(x = Day, y = Bs)) +
theme_bw()
## `summarise()` has grouped output by 'Description'. You can override using the
## `.groups` argument.
Myxococcota - вроде бы тоже, но как нам известно могут быть целлулотитиками
ps.m %>%
filter(Phylum == "Myxococcota") %>%
group_by(Description, Day) %>%
summarise(Bs = sum(Abundance)) %>%
ggplot() +
geom_boxplot(aes(x = Day, y = Bs)) +
theme_bw()
## `summarise()` has grouped output by 'Description'. You can override using the
## `.groups` argument.
Археи появляются тоже на поздних стадиях - вообще я бы хотел бы опять
развить тему важности азотного метаболизма на поздних стадиях
разложения.
Оч хочется метагеном, но не этот.
Кроме того хочется отметить, что эти минорные группы возникают на
D12
ps.m %>%
filter(Phylum == "Crenarchaeota") %>%
group_by(Description, Day) %>%
summarise(Bs = sum(Abundance)) %>%
ggplot() +
geom_boxplot(aes(x = Day, y = Bs)) +
theme_bw()
## `summarise()` has grouped output by 'Description'. You can override using the
## `.groups` argument.
Gammaproteobacteria - они треть кластера blue
ps.m %>%
filter(Class == "Gammaproteobacteria") %>%
group_by(Description, Day) %>%
summarise(Bs = sum(Abundance)) %>%
ggplot() +
geom_boxplot(aes(x = Day, y = Bs)) +
theme_bw()
## `summarise()` has grouped output by 'Description'. You can override using the
## `.groups` argument.
Все филы - представлененось логорифмирована по основанию 2
Отдельные точки - суммы абсолютных значений ридов по дням
Логорифмированы уже суммы, а не отдельные филотипы
#select only major phylums
top_phylum <- ps.m %>%
count(Phylum) %>%
arrange(desc(n)) %>%
top_n(10) %>%
pull(Phylum)
## Selecting by n
ps.m %>%
filter(Phylum %in% top_phylum) %>%
mutate(
Phylum = as.character(Phylum),
Class = as.character(Class),
phylum = ifelse(Phylum == "Proteobacteria", Class, Phylum)
) %>%
group_by(Description, Day, phylum) %>%
filter(!is.na(phylum)) %>%
summarise(Bs = log2(sum(Abundance))) %>%
ggplot(aes(x = Day, y = Bs)) +
geom_boxplot(fill="#4DBBD5B2", alpha=0.4) +
theme_bw() +
facet_wrap(~ phylum)
## `summarise()` has grouped output by 'Description', 'Day'. You can override
## using the `.groups` argument.
## Warning: Removed 46 rows containing non-finite values (stat_boxplot).